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the  FY  77  HPR  Work  Program  on  this  Study.   A  controlled  laboratory 
study  will  be  proposed  for  verification  and  fine-tuning  of  the 
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HIGHLIGHT  SUMMARY 


An  analytical  tool,  based  on  the  finite  element  method,  has  been 
developed  to  analyze  buried  culvert  problems  in  a  realistic  fashion. 

Segments  of  a  curved  bar  with  three  degrees  of  freedom  (normal, 
tangential  and  rotational)  at  each  end  have  been  used  to  simulate  a 
thin  pipe  where  nodal  moments  are  important.   Triangular,  isopara- 
metric elements  with  one  curved  boundary  (to  fit  the  shape  of  pipe), 
and  three  midside  nodes  have  been  used  to  represent  the  soil.   A 
special  type  of  'interaction'  element  with  zero  thickness  has  been 
used  between  pipe  and  soil  to  simulate  interface  behavior,  including 
slip. 

Nonlinear,  anisotropic  soil  properties  have  been  accounted  for. 
Actual  test  data  are  used  as  input  for  soil  properties;  the  data  are 
interpolated  using  a  cubic  spline  function,  incremental  Young's 
modulus  and  Poisson's  ratio  are  calculated  in  terms  of  octahedral 
normal  and  shear  stresses  and  the  result  is  stored  in  the  computer. 
'Interaction'  elements  offer  high  resistance  to  movements  in  the 
normal  and  tangential  directions  until  the  stress  ratio  (shear  to 
normal  stress)  exceeds  a  limiting  value.   Once  the  ratio  is  exceeded, 
stiffness  in  the  shear  direction  is  reduced  to  zero.   Actual  construction 
sequences  can  be  simulated  and  'no-tension'  analysis  is  incorporated  if 
the  soil  and/or  the  soil-pipe  interface  cannot  resist  tension. 
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A  computer  program  was  written  to  include  all  the  aspects 
mentioned  above,  and  some  example  problems  were  solved  to  demonstrate 
its  versatility  and  to  investigate  the  influence  of  such  factors  as 
non-linear  soil  properties,  relative  stiffness  of  pipe  and  soil, 
inclusion  of  weak  materials  near  the  spring  line,  and  construction 
procedures.   It  was  found  that  for  circular  corrugated  metal  pipe 
buried  in  granular  soil: 

(a)  maximum  circumferential  thrust  in  the  pipe  depends  mainly  on  pipe 
diameter  and  height  of  fill,  although  interface  friction  has  some 
influence. 

(b)  the  effect  of  nonlinear  soil  properties  on  culvert  performance 
must  be  accounted  for  but,  for  high  fills  placed  and  loaded 
symmetrically,  neither  construction  sequence  nor  partial  slip  at  the 
soil-pipe  interface  had  significant  effects. 

(c)  Marston-Spangler  soil  modulus  E'  can  not  be  used  as  a  soil 
parameter  in  rational  predictions  of  culvert  performance. 

The  study  did  not  consider  time-dependent  soil  properties, 
prestress  of  soil  during  compaction,  or  cracking  or  yielding  of  the 
pipe,  although  the  latter  can  readily  be  accounted  for.   The 
present  program  is  applicable  only  for  analysis  of  two-dimensional 
problems  transverse  to  the  pipe  in  which  the  state  of  stress  in  the 
soil  mass  does  not  approach  failure. 


CHAPTER  I 
INTRODUCTION 

Culverts  and  conduits  under  highway  and  railway  em- 
bankments are  commonly  used  structures.   Depending  upon  the 
purpose,  they  may  be  large  or  small,  deep  or  shallow,  rigid 
or  flexible.   Culverts  can  be  classified  in  two  broad  cat- 
egories, based  on  their  relative  stiffness:   (1)   rigid  and 
(2)   flexible.   Rigid  culverts  are  those  for  which  the 
diametral  change  under  load  prior  to  rupture  is  assumed  to 
be  too  small  to  influence  the  resulting  distribution  of  soil 
pressure.   Before  cracking,  concrete  pipe  culverts  are  usu- 
ally considered  to  be  rigid.   Flexible  culverts  are  designed 
on  the  basis  that  sufficient  deflection  of  the  culvert  will 
occur  to  mobilize  additional  lateral  resistance  from  the 
surrounding  soil  mass.   Accordingly,  for  rigid  culverts,  the 
"Marston  Theory"  (Marston  and  Anderson,  1913;  Marston,  19  30) 
is  commonly  used  for  design  while  flexible  culverts  are  de- 
signed using  the  well-known  "Iowa  Formula",  (Spangler,  1962; 
FHWA,  1970) .   For  rigid  culverts  bedded  within  an  earth  fill, 
the  Marston  approach  assumes  that  the  total  load  applied  to 
the  horizontal  projection  of  a  conduit  is  equal  to  the 
weight  of  the  fill  above  the  conduit,  modified  by  a  factor 


which  purports  to  account  for  the  relative  displacement  of 
the  prism  of  soil  above  the  conduit  with  respect  to  the  sur- 
rounding soil.   For  design  purposes,  Spangler  (1962)  recom- 
mends that  this  factor  be  equal  to  or  greater  than  unity,  de- 
pending upon  the  foundation  condition.   The  rigid  culvert  is 
designed  to  resist  this  uniform  vertical  load  in  addition  to 
those  transmitted  by  external  dead  and  live  loads.   But  in 
determining  the  lateral  pressure,  which  assists  in  supporting 
the  culvert,  no  account  is  taken  of  the  stresses  in  the  soil 
surrounding  the  pipe.   For  flexible  culverts,  the  design  is 
based  on  limiting  the  change  in  vertical  diameter  of  the 
pipe,  as  well  as  providing  adequate  resistance  against  crit- 
ical buckling  of  the  pipe  wall,  or  failure  of  the  longitudi- 
nal (riveted,  welded  or  bolted)  seams  (White  and  Layer, 
1960;  Spangler,  1962;  FHWA,  1970)  .   The  lateral  resistance 
offered  by  the  soil  is  characterized  by  a  so-called  modulus 
of  soil  reaction  (E1),  which  is  defined  (Spangler,  1962)  as 
the  ratio  of  the  lateral  stress  to  the  displacement  at  the 
spring  line  of  the  pipe,  times  the  pipe  radius.   Values  of  E' 
are  recommended  based  on  field  experience  and  model  tests. 

The  semi-empirical  nature  of  present  design  methods, 
expecially  in  evaluating  the  soil-structure  interaction, 
renders  them  difficult  to  extrapolate  to  situations  other 
than  those  which  have  been  considered  in  the  past.   This  fea- 
ture is  of  particular  importance  in  view  of  the  recent  trend 
toward  high  fills  (Davis  and  Bacher,  1968) ,  and  large 


diameter  culverts  under  shallow  fills.   Extrapolation  be- 
comes even  more  difficult  when  new  materials  such  as  plas- 
tics and  composites  are  introduced. 

The  main  purpose  of  this  study  is  to  develop  an 
analytical  design  tool  which  can  be  used  to  simulate  actual 
conditions  realistically,  and  to  predict  the  performance  of 
culverts  embedded  in  soil  under  different  construction  and 
loading  situations.   All  efforts  have  been  directed  towards 
representative  simulation  of  pipe,  soil,  and  soil-pipe 
interaction.   For  this  purpose  the  technique  of  finite  ele- 
ment analysis  was  chosen.   Flexible  pipes  were  simulated  by 
segments  of  a  curved  bar.   Isoparametric  triangular  elements 
with  stress  dependent,  nonlinear  soil  properties  were  used 
to  represent  the  soil  medium.   A  special  type  of  interaction 
element  was  introduced  between  soil  and  pipe  which  permits 
slip  to  occur  between  soil  and  pipe  at  shear  stresses  less 
than  the  strength  of  the  soil.   Incremental  analysis  to 
simulate  construction  in  steps  has  been  accomodated. 

A  finite  element  computer  program  has  been  written 
which  incorporates  all  the  factors  mentioned  above.   This 
program  has  been  used  in  a  few  example  problems  to  show  its 
capability  for  the  purpose  of  predicting  culvert  perform- 
ance in  practical  situations.   Examples  were  selected  to  show 
the  effects  of  various  soil  parameters  and  different  types  of 
loadings. 


The  finite  element  technique  developed  in  this  study 
was  used  to  investigate  non-linear  vs.  linear  soil  properties, 
the  effects  of  construction  in  layers,  inclusion  of  stiff er 
or  softer  material  at  different  locations  on  the  pipe,  slip 
between  pipe  and  soil  and  various  combinations  of  super- 
imposed loadings.   It  is  believed  that  the  significant  as- 
pects of  culvert  behavior  have  been  realistically  repre- 
sented and  that  it  should  be  possible  to  predict  performance 
under  a  wide  variety  of  field  conditions.   Comparison  of  such 
predictions  with  measured  performance  (from  controlled  lab- 
oratory tests  and  field  installations)  are  needed  to  estab- 
lish the  practical  utility  of  the  newly  developed  method  to 
predict  actual  field  performance. 


CHAPTER  II 
LITERATURE  REVIEW 

2. 1   Current  Design  Practice 

In  this  section  current  design  approaches  of  differ- 
ent types  of  culverts  are  reviewed  with  emphasis  on  the 
assumptions  made  and  their  validity. 

Concrete  Culverts:   (Rigid  culverts)  -  California  Cul- 
vert Practice  (1944)  mentions  the  importance  of  differential 
settlement  between  the  pipe  and  surrounding  soil  and  its 
effect  on  the  loads  carried  by  the  culvert.   For  a  rigid  cul- 
vert on  unyielding  foundation,  the  earth  alongside  the  cul- 
vert moves  downwards  relative  to  the  material  in  the  prism 
over  the  culvert.   This  causes  load  to  be  transferred  to  the 
prism  over  the  culvert  which  will  result  in  an  excess  load 
over  the  culvert.   When  the  culvert  settles  or  deflects  an 
amount  equal  to  the  settlement  of  the  surrounding  material, 
the  load  on  culvert  will  be  the  same  as  the  weight  of  the 
prism  of  material  directly  above  it.   This  mechanistic  pic- 
ture of  behavior  is  very  realistic  but  the  amount  of  relative 
settlement  and  their  effects  must  be  evaluated  by  qualitative 
judgment.   Once  this  is  done  certain  specific  design 
approaches  are  available.   For  example,  Spangler,  (1960) 


presents  a  formulation  for  determining  vertical  loads  for  un- 
yielding culverts,  which  he  named  'projecting'  culverts. 
Townsend  (1963)  presented  design  charts  for  estimating  loads 
on  pipes  for  various  projecting  conditions.   His  specifica- 
tions estimating  loads  of  ASTM  standard  (ASTM,  C76-59T)  pipes 
are  based  mainly  on  Marston's  formula,  and  his  charts  merely 
transform  Marston's  formula  into  graphical  form. 

Marston  Theory  -  Marston's  theory  is  often  used  to  de- 
termine vertical  loads  on  flexible  or  rigid  culverts.   The 
mechanism  and  assumptions  of  this  theory  is  illustrated  in 
Figure  (2.1),  in  which 

P   =  External  vertical  load 

V  =  Vertical  force  on  a  horizontal  plane 

dV  =  Increment  of  vertical  force 

B   =  External  diameter  of  pipe 
c 

B-,  =  Horizontal  width  of  ditch 
d 

H   =  Height  of  fill  on  top  of  pipe 
h  =  Variable  distance  from  top  of  soil 
dh  =  Increment  of  h 

K   =  Coefficient  of  active  earth  pressure 
u'  =  tan  <j> '  =  Coefficient  of  friction  between  fill 
material  and  side  of  ditch 

Y  =  Unit  weight  of  fill. 

For  equilibrium  of  vertical  forces  in  the  element  of  thick- 
ness dh,  the  following  equation  can  be  written  per  unit 
length  of  pipe  (assuming  P=0) : 


pwr 


dh 


H 


V 
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FIGURE   2.1     FREE    BODY    DIAGRAM     FOR    DITCH 
CONDUIT.  (After   Spongier  ,  I960) 


V+dV+2Ky'^-  dh  =  V+yB,  dh  (2.1) 

Bd  d 


Solving  for  V, 

-2Ky'  (£-) 

V  =  YBd  Tfi "  Cd  ^Bd  (2'2) 

On  substitution  of  h  =  H,  total  vertical  load  on  the 
pipe  can  be  evaluated.   According  to  Mars ton  theory  equation 
(2.2)  represents  the  total  vertical  load  on  pipe  for  rigid 
culverts  with  relatively  compressible  side  fills.   For  flex- 
ible culvert,  with  thoroughly  compacted  side  fills  which  have 
essentially  the  same  degree  of  stiffness  as  the  pipe,  V  is 
modified  as: 

V  =  C ,  y  B   B,  (2.3) 

d    c   d 

Spangler  (1963)  points  out  that  in  actual  cases  V 
will  lie  between  the  values  given  by  equations  (2.2)  and 
(2.3).   This  formulation  is  logical,  but  leaves  some  ques- 
tions unanswered. 

(1)  For  how  wide  a  ditch  are  the  formulas  applica- 

ble?   Spangler  (1963)  presented  charts  for  C,  values  for  — 

a  ud 

values  ranging  up  to  15,  which  shows  values  of  C,  increases 

for  —  up  to  about  8.   He  mentioned  the  width  of  the  ditch  as 

Bd 
"relatively  narrow"  and  he  did  not  present  any  criteria  for 

establishing  limiting  values  of  B^. 

(2)  Is  it  proper  to  consider  active  earth  pressure 
condition  with  the  wall  of  trench  without  knowing  the 
deflection? 


(3)  What  will  be  the  changes  in  the  load  on  pipe  for 
different  specifications  for  degree  of  compaction  of  fill 
surrounding  the  pipe? 

(4)  Will  there  be  any  difference  in  load  if  there  is 
slip  between  pipe  and  soil? 

(5)  Will  the  load  P  be  distributed  to  the  pipe  ac- 
cording to  the  Boussinesq  solution,  as  proposed  by  Marston? 

These  and  other  factors  were  accounted  for  by 
Spangler  (1963)  by  defining  specific  restraint  conditions 
(positive  projection,  negative  projection,  imperfect  ditch, 
settlement  ratio)  and  preparing  design  charts  based  on  model 
tests  and  field  experiences. 

After  the  load  on  the  pipe  has  been  computed,  selec- 
tion of  pipe  is  based  on  either  of  two  criteria,  (i)  0.01 
inch  crack  strength, or  (ii)  ultimate  pipe  strength   divided 
by  an  appropriate  value  of  Factor  of  Safety.   Recommended 
values  of  Factory  of  Safety  are  1.5  for  the  first  criteria 
and  1.33  for  the  second.   Pipe  strengths  are  determined  by 
conventional  three-edge-bearing  tests. 

Concrete  pipe  culverts  do  not  collapse  when  cracking 
is  initiated,  and  as  further  cracking  develops   the  rigidity 
of  the  pipe  changes  sufficiently  to  redistribute  the  soil 
pressures  on  the  pipe.   Important  research  to  predict  per- 
formance of  cracked  concrete  pipe  is  underway  at  Northwestern 
University,  (Krizek  et.  al.  (1971) ) . 
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Flexible  Culverts  -  Wolf  and  Townsend  (1970)  proposed 
that  the  following  criteria  should  guide  the  design  of  flex- 
ible culverts. 

(1)  deflection  or  flattening  of  pipe 

(2)  buckling  of  pipe  wall 

(3)  longitudinal  seam  strength 

(4)  handling  and  installation  strength. 
Deflection  criterion  -  Change  in  nominal  pipe  diameter 

for  flexible  culverts  can  be  determined  using  the  "Iowa  Form- 
ula" (Spangler,  1960)  : 

DTKW   R 
Ax  =  _i_E (2.4) 

EI  +  0.061  RE1 


where 


Ax  =  horizontal  deflection  of  pipe  in  inches, 

assumed  equal  to  vertical  deflection 

D  =  deflection  lag  factor,  (1.25  for  E'  =  1400  psi 
L 

or  more   1.5  for  E'  =  700  psi) 
R  =  radius  of  pipe  in  inches 
K  =  bedding  constant,  depending  on  the  bedding  angle 

cc(0.11  for  a  =  0  to  0.083  for  a  =  90°) 

W  =  vertical  load  on  pipe  as  predicted  by  Marston 
c 

theory 
E  &  I  =  Young's  modulus  and  moment  of  inertia  of  the 

pipe 
E'  =  modulus  of  "soil  reaction",   in  pounds 
per  square  inch 
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This   formula  has  been  derived  based  on  the  theory  of 
curved  beams  and  arches  on  elastic  foundation  and  assumes  a 
parabolic  distribution  of  horizontal  stresses  at  side  of  pipe. 
The  Iowa  Formula  is  based  on  an  assumed  distribution  of  load 
around  pipe,  and  parameters   DT ,   K   and  E*  are  empirically 
derived.   Parmelee  and  Corotis  (1972)  took  a  closer  look 
into  the  assumptions  and  validity  of  Iowa  Formula.   They 
studied  data  from  18  tests  and  conducted  statistical  analysis 
of  the  three  important  parameters   D  ,   K   and  E'. 
They  concluded  that  there  is  no  rational  basis  for 
assigning  any  value  for  E1  which  is  greater  than  the  modulus 
value   (E)   of  the  medium.   They  also  pointed  out  that  ex- 
treme caution  should  be  exercised  in  choosing  values  for  the 
parameters  in  the  Iowa  Formula,  including  values  of  W  . 


The  deflection  criterion  arbitrarily  specifies  fail- 
ure at  5  percent  deflection  of  nominal  pipe  diameter.   Ex- 
periments suggest  that  a  20  percent  change  in  vertical  diameter 
is  necessary  to  develop  plastic  hinges  in  corregated  metal 
pipes,  provided  local  buckling  does  not  occur.   In  such  cases, 
the  deflection  criterion  has  a  nominal  Factor  of  Safety  of  4. 

Buckling  criterion  -  This  criterion  provides  for  the 
design  of  pipe  based  on  the  wall  thickness  required  for  a 
limiting  buckling  stress  which  takes  into  account  the  re- 
straining effect  of  soil-structure  interaction  around  pipe. 
This  restraining  effect  is  solely  due  to  confinement  of  the 
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pipe  in  soil.   If  the  pipe  cross-section  is  considered  as 
ring,  then  the  ring  compression  becomes  critical  for  buck- 
ling when, [Wolf  and  Townsend  (1970).] 

f  2 
F   >  f   -  JL.  (]LA)2  (2.5) 

b  -  u    48E  l  r  ' 


for  diameter  less  than  =-  /  24E 


/ 


f 
u 


where 


F,  =  buckling  stress,  psi 


f   =  minimum  tensile  strength,  psi 

k   =  soil  stiffness  coefficient  (0.22  for  E1  =  1400 
psi  or  above  and  0.44  when  E'  =  700  psi.) 

d  =  nominal  diameter  of  pipe,  inch 

r   =  radius  of  gyration  of  pipe  wall,  inch. 

E   =  Young's  modulus  of  pipe  material,  psi 

Design  for  buckling  is  accomplished  by  limiting  the 
ring  compression  thrust  to  the  buckling  stress  multiplied  by 
the  conduit  wall  area,  divided  by  Factor  of  Safety. 

Longitudinal  seam  strength  and  handling  stresses  are 
not  considered  in  this  study,  although  their  relevancy  to  the 
total  problem  is  recognized.   It  might  be  mentioned,  however, 
that  yielding  (as  opposed  to  rupture)  of  a  seam  is  often 
beneficial  as  it  is  equivalent  to  an  increase  in  pipe  flexi- 
bility without  necessarily  reducing  the  buckling  resistance. 

It  is  seen  that  the  Iowa  Formula  assumes  both  soil 
and  pipe  as  linear,  elastic  and  homogeneous  materials  during 
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the  entire  loading  range,  hence  values  of  Wn,  DT ,  K  and  E1 
are  difficult  to  establish  for  wide  ranges  in  design  condi- 
tions.  Several  manuals  provide  numerical  values  of  these 
parameters  based  mostly  on  past  experience.   Culvert  design 
methods  suggested  by  pipe  manufacturers  (e.g.,  Republic 
Steel,  ARMCO)  follow  simplified  procedures  as  follows: 

(1)  Dead  load  -  the  total  weight  of  the  prism 
of  soil  directly  above  the  pipe. 

(2)  Live  load  -  specified  by  empirically  derived 
charts. 

(3)  Ring  thrust  -  vertical  pressure  multiplied  by 
the  pipe  radius. 

From  the  value  of  ring  thrust  the  required  sectional  proper- 
ties of  pipe  are  chosen.   None  of  these  manuals  attempt  to 
evaluate  soil-structure  interaction,  nonlinear  soil  proper- 
ties, or  installation  sequence,  but  most  mention  the  impor- 
tance of  compacting  the  backfill  materials  close  to  pipe 
wall. 

2. 2   Literature  Review 

Krizek  and  Parmelee  (1968)  prepared  an  extensive 
bibliography  on  analysis,  design  and  installation  of  highway 
culverts,  which  covers  publications  from  year  1900  to  1968. 
A  review  of  selected  recent  studies  follows. 

Davis  (1966)  reported  results  of  a  field  study  of  a 
large  reinforced  concrete  arch  culvert  buried  under  an 
embankment  200  feet  high  in  natural  soil.   This  report 
concluded  that 
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(1)  stresses  at  any  depth  are  proportional  to  the 
fill  height 

(2)  soil  stresses  remain  essentially  constant  after 
the   fill  is  completed 

(3)  side  walls  of  arch  moved  inwards  while  the  crown 
moved  upwards  during  back  filling 

(4)  a  layer  of  organic  material  around  the  pipe  re- 
duced vertical  pressure  to  about  half  of  that 
expected 

(5)  distribution  of  soil  pressure  around  culvert  was 
vastly  different  from  the  values  predicted  from 
a  model  prepared  by  Brown  (1967). 

Brown,  Green  and  Pawsey  (1968)  conducted  tests  on 
culverts  under  high  fill  with  soft  inclusions  (straw)  at 
different  locations  around  pipe,  which  showed  significant 
difference  in  stress  distribution  around  pipe.   Theoretical 
results  using  incremental  finite  element  analysis  did  not 
predict  the  nature  of  the  experimental  results. 

Watkins  and  Moser  (1971)  presented  simple  design 
criteria  for  flexible  pipes  based  on  test  data  from  a  number 
of  field  and  laboratory  studies.   They  assumed  a  section  of 
a  pipe  as  a  ring  and  related  the  ring  compression  strength, 
f  with  the  applied  stress  by  the  simple  formula. 

PD  -   c  (0    M 

2A  "  W  (2'6) 
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where 

P  =  vertical  soil  pressure  determined  by  adding 

pressure  on  pipe  due  to  weight  of  prism  of  soil 
and  that  produced  by  any  imposed  live  load. 

D  =  nominal  diameter  of  pipe 

A  =  cross-sectional  area  of  pipe 

N  =  Factory  of  Safety 
Equation  (2.6)  states  that  for  limiting  equilibrium,  applied 
stress  is  equal  to  the  strength  of  pipe  divided  by  adequate 
Factor  of  Safety.   Based  on  this,  and  the  results  of  labora- 
tory tests,  the  authors  presented  a  design  chart,  which  in- 
cludes the  effect  of  degree  of  compaction  of  the  surrounding 
soil.   This  procedure  does  not  take  into  account  soil- 
structure  interaction  due  to  arching  and  terminal  shear  be- 
tween pipe  and  soil.   Spangler  (1971)  points  out  that  the 
model  proposed  by  Watkins  and  Moser  (1971)  ignores  the 
moments  developed  in  the  pipe  which,  when  neglected,  is  not 
representative  of  pipe  behavior.   Spangler  (1971)  concluded 
that  the  similarity  between  the  field  test  results  and 
equation  (2.6)  is  coincidental.   However,  Neilson  (1972) 
pointed  out  that  the  test  results  by  Watkins  and  Moser  (1971) 
have  good  agreement  with  the  results  given  by  elastic  theory 
using  constrained  soil  modulus  values,  and  proposed  that 
agreement  of  laboratory  and  field  results  by  Watkins  and 
Moser  (1971)  may  be  due  to  same  degree  of  confinement  in  both 
cases. 
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Allgood  and  Takahashi  (1972)  have  shown  the  similar- 
ity between  Iowa  Formula  and  the  solution  for  an  elastic 
ring  in  an  elastic  medium,  and  presented  a  semi-empirical 
procedure  for  culvert  design  whose  principal  advantages  are: 

(1)  Arching  of  soil  can  be  considered 

(2)  Effects  of  backpacking  can  be  considered 

(3)  Moment  in  pipe  can  be  taken  into  account 

(4)  Several  modes  of  failure  are  taken  into  account. 
This  procedure  assumes  both  soil  and  pipe  as 

linearly  elastic  material  and  ignores  potential  slip  at  the 
soil-structure  interface. 

Allgood  and  Takahashi  (1972)  also  presented  results 
of  three  dimensional  finite  element  analysis  of  culvert 
problems.   Changes  in  elastic  modulus  were  considered  by 
specifying  values  at  different   depths.   Poisson's  ratio  was 
kept  unchanged  at  a  given  value.   Possibility  of  slip  between 
pipe  and  soil,  non-linear  soil  properties,  sequence  of  con- 
struction, and  variations  in  state  of  compaction  around  pipe 
were  neglected.   As  a  result,  the  authors  acknowledged  that 
large  differences  can  be  expected  between  their  predictions 
and  field  performance. 

Kirkland  and  Walker  (1972)  performed  model  studies  of 
culverts  on  6  inch  diameter  steel  pipes  with  different  wall 
thickness  so  as  to  consider  both  rigid  and  flexible  culverts. 
Soil  properties  were  modeled  in  two  ways:   (1)   use  of  bulk 
and  shear  modulus  updated  for  state  of  volumetric  and  shear 
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strains,  (2)   by  fitting  second  degree  curves  for  bulk  and 
shear  modulus.   Behavior  for  soil-structure  interaction  was 
simulated  by  two  curves,  one  for  initial  near-straight  line 
portion  and  the  second  for  the  portion  where  the  curve  be- 
gins to  flatten.   Comparison  of  predicted  and  experimental 
results  always  show  significant  deviation.   This  may  be 
attributed  to  the  following  sources  of  error   (1)   assumption 
of  isotropic  behavior  of  soil,  (2)   linear  displacement  func- 
tion for  finite  elements,  (3)  homogeneous  soil  properties. 
The  authors  indicated  the  need  for  more  research  on  the  effects 
of  slip  at  the  soil-pipe  interface. 

Howard  and  Metzger  (1973) ,  Howard  (1973)  conducted  a 
large  number  of  laboratory  and  field  tests  on  a  variety  of 
pipes  which  includes  soil,  fiberglass,  reinforced  resin  and 
reinforced  plastic  mortar  for  pipe  materials.   Results  show 
significant  effects  of  various  types  of  bedding  and  compac- 
tion around  pipe.   For  well  compacted  backfills,  pipes  of 
all  variety  of  materials  show  similar  response.   But  for 
poorly  compacted  backfills,  the  results  vary  over  a  wide 
range.   The  spread  of  results  can  be  attributed  to  the  vari- 
ation of  friction  between  soil  and  different  pipe  materials, 
which  has  significant  effects  for  relatively  low  soil  density. 
It  also  shows  that  behavior  of  pipes  with  reinforced  lam- 
inates significantly  differs  from  that  of  steel  pipes.   These 
aspects  need  further  study.   Davis,  Bacher  and  Obermuller 
(1974)  conducted  field  tests  of  an  under-designed  reinforced 
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concrete  pipe  culvert  of  large  diameter  (84")  in  a  high 
embankment  (136  ft) .   Main  purpose  of  this  study  was  to  in- 
vestigate the  structural  behavior  of  several  methods  of 
backfilling  and  bedding.   It  was  found  that  with  proper 
attention  to  bedding  and  backfilling  procedures,  in  partic- 
ular, providing  adequate  lateral  support  by  placing  the  pipe 
in  a  trench  and  surrounding  it  by  a  well-compacted  high 
quality  structural  backfill,  culverts  of  much  lesser  thickness 
than  normal  can  be  used  satisfactorily .  In  a  related  paper, 
Davis  and  Bacher  (1974)  concluded  that  for  embankments  of  appre- 
ciable height,  earth  pressure  loads  are  of  far  more  consequence 
than  any  external  live  loads.   The  authors  expressed  their 
concern  over  the  design  and  backfill  procedures  specified 
by  California  Division  of  Highways  Specifications.   They  in- 
dicated a  need  for  future  research  in  the  following  areas: 

(1)  Soil-pipe  interface  around  pipe 

(2)  Different  types  and  thickness  of  bedding 

(3)  Use  of  a  layer  of  polystyrene  surrounding  the 
pipe  to  reduce  friction 

(4)  Revision  of  current  tables  and  design  charts. 
Dar  and  Bates  (1974)  presented  a  method  for  analysis 

of  culvert  problems  based  on  theory  of  elasticity  for  ideal 
linear  elastic,  isotropic  and  homogeneous  soil  and  pipe 
properties.   Their  derivation  treats  the  pipe  as  a  circular 
elastic  material  and  the  medium  as  an  elastic  plate  with  a 
circular  hole.   Stresses  in  the  pipe  are  determined  from  a 
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harmonic  displacement  function  and  the  stresses  in  the  medium 
is  determined  from  traction  stresses  around  the  opening. 
Matching  the  stresses  and  displacements  at  the  boundary  be- 
tween pipe  and  soil,  and  permitting  no  slip  at  the  boundary, 
a  solution  to  the  problem  is  obtained.   Experimental  results 
for  concrete  pipes  were  compared  with  those  predicted  by  the 
formula.   Results  show  that  the  method  provides  an  upper 
bound  for  radial  stresses  in  thin  pipes,  and  a  lower  bound  for 
thick  pipes.   This  distinct  difference  in  behavior  may  be 
due  to  the  interface  conditions  and  simplifying  assumptions 
used  in  the  analysis.   The  approach  is  of  value  in  estimat- 
ing upper  and  lower  bounds  in  the  pipe  stresses  but  more 
comparisons  with  field  measurements  are  needed. 

A  long-range  research  project  concerning  soil- 
structure  interaction  of  concrete  culverts  is  underway  at 
Northwestern  University,  Evanston,  Illinois.   This  research 
program  consists  of  two  parts,  (i)   a  series  of  full  scale 
field  installation  of  concrete  culverts  including  instrumen- 
tation, (ii)   development  of  finite  element  program  to  pre- 
dict behavior.   Nonlinear  properties  for  both  soil  and  pipe 
have  been  considered.   Cracking  of  concrete  pipes  will  be 
modelled  in  this  research  program.   Parmelee  (1973)  mentioned 
that  decisions  about  the  soil-structure  interface  will  be  made 
after  developing  considerable  data  about  their  interaction 
through  a  series  of  laboratory  and  field  tests. 
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2 . 3   Need  for  Further  Research 

Both  flexible  and  rigid  culverts  are  commonly  in  use. 
Their  size  varies  from  several  inches  to  tens  of  feet 
(Fisher,  1969).   Shapes  vary  over  a  wide  range  which  includes 
circular,  elliptical,  bullet-shape,  arch,  and  combinations  of 
different  shapes.   Depth  of  cover  vary  from  several  feet  to 
several  hundred  feet.   Methods  of  construction  include 
trenching,  bedding,  backfilling  and  compaction.   Different 
backfill  materials,  ranging  from  loose  straw   to  foamed 
plastics  are  used  around  the  pipe  to  change  stress  distri- 
bution patterns.   The  soil  medium  varies  from  coarse  granu- 
lar material  to  plastic  clay.   Materials  used  for  the  pipe 
include  reinforced  concrete,  prestressed  concrete,  plate 
steel,  corrugated  steel,  aluminum, clay  tile,  reinforced 
plastic,  etc.  Superimposed  loads  may  vary  over  a  wide  range. 
Thus,  the  number  of  variables  encountered  in  a  real  culvert 
problems  is  enormous,  and  can  only  increase  in  the  future. 
This  imposes  severe  restrictions  on  approaches  that  rely  on 
empirically  derived  parameters  and  simplified  design  proce- 
dures. 

It  is  proposed  in  this  study  to  develop  an  analytical 
tool  which  will  permit  a  realistic  study  of  culvert  problems 
with  a  minimum  number  of  restrictive,  simplifying  assumptions, 
This  tool  is  intended  to  be  capable  of  handling  generalized 
soil  properties,  soil-culvert  interaction  including  slip  at 
the  interface,  and  any  combination  of  construction  conditions 
and  superposed   loadings. 
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CHAPTER  III 
FINITE  ELEMENT  METHOD 

3. 1   Historical  Background. 

The  finite  element  method  of  analysis  is  a  well- 
known  tool  in  many  branches  of  engineering  science  and  tech- 
nology.  Many  complex  problems  can  be  solved  satisfactorily 
by  the  use  of  this  method. 

The  concept  of  Ritz  (Ritz,  1909)  and  Raleigh's 
(Strutt  ,  1870)  method  of  variational  analysis  has  been  in 
existence  for  long  time.   In  the  Rayleigh-Ritz  method,  the 
expression  for  total  potential  energy  of  a  system  is  formu- 
lated and  the  displacement  pattern  is  assumed  to  vary  with  a 
finite  number  of  undetermined  parameters.   A  set  of  simultan- 
eous equations  are  formulated  minimizing  the  total  potential 
energy  with  respect  to  the  displacement  parameters.   So  far 
both  the  finite  element  and  Rayleigh-Ritz  methods  are 
identical.   In  the  Ritz  procedure,  the  displacement  pattern 
is  specified  in  the  whole  region,  whereas  in  the  finite  element 
method  the  displacements  (or  loads)  are  specified  at  nodal 
points  of  discrete  elements. 

The  organized  matrix  approach  to  structural  analysis 
based  upon  the  finite  element  idealization  was  prescribed 
clearly  by  Turner,  et.  al.  (1956)  in  their  pioneering  paper. 
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In  this  paper  important  concepts  of  compatibility  at  the 
boundary  between  two  elements  in  terms  of  nodal  displacements, 
nodal  equilibrium,  and  determination  of  stiffness  influence  co- 
efficients were  presented  in  step  by  step  development.   The 
main  criteria  of  variational  methods  in  linear  static  struc- 
tural mechanics  is  the  principle  of  minimum  potential  energy 
(Washizu,  1968)  ,  which  provides  a  basis  for  the  direct  formu- 
lation of  element  stiffness  equations.   The  paper  by  Turner, 
et.  al.  did  not  make  reference  to  variational  considerations, 
rather,  it  took  a  "direct"  approach,  in  which  direct  consid- 
eration was  given  to  conditions  of  equilibrium  and  compati- 
bility.  Consequently,  the  paper  stimulated  numerous  efforts 
to  develop  formulations  for  stiffness  influence  coefficients 
for  a  variety  of  different  types  of  elements  (Turner,  1959; 
Clough,  1960;  Melosh,  1961;  Melosh,  1963;  Clough  and  Tocher, 
1965)  . 

In  the  last  decade,  research  in  finite  element  meth- 
ods were  directed  mainly  to  new  techniques  of  analysis  for 
complex  structures  with  nonlinear  material  properties,  heat- 
flow  and  seepage  (with  moving  boundaries) ,  large  displace- 
ments, etc.  (Mallett  and  Marcal,  1968;  Przemieniecki ,  1968; 
Zienkiewicz,  1971;  Desai  and  Abel,  1972).   The  method  has 
been  extended  to  deal  with  three-dimensional  geometries  using 
sub  and  superparametric  elements,  curvilinear  elements,  and 
account  has  been  taken  of  inertia  effects. 
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3. 2   Requirements  for  Finite  Element  Methods 

As  discussed  in  Chapter  II,  the  behavioral  aspects 
of  culverts  can  be  analyzed  with  confidence  provided  each 
component  of  the  structure  such  as  soil,  pipe,  and  soil- 
structure   interface   can  be  modeled  realistically  using  a 
minimum  number  of  restrictive  assumptions.   To  achieve  this, 
the  first  step  is  to  idealize  all  components  of  the  problem  by 
an  appropriate  geometric  shape.   A  good  choice  of  element 
shape  should  consider  the  following  factors. 

(1)  The  elements  can  be  assembled  to  form  the  com- 
plete structure  with  correct  shapes  at  the  boundaries. 

(2)  The  element  should  have  a  sufficient  number  of 
degrees  of  freedom  to  simulate  the  actual  structural  response, 

(3)  The  element  should  not  be  ill-conditioned,  i.e., 
convergence  is  too  slow  or  is  prevented. 

Once  the  element  shape  is  chosen,  the  next  step  is 
to  choose  the  displacement  function  (also  known  as  shape 
function) .   Careful  consideration  of  the  following  factors 
are  essential  for  the  choice  of  shape  function. 

(1)  The  displacement  function  should  ensure  minimum 
potential  energy. 

(2)  Continuity  is  maintained  at  the  boundary  between 
adjacent  elements. 

(3)  Straining  of  an  element  due  to  rigid  body  motion 

is  precluded. 

(4)  Criteria  for  convergence  to  the  true  solution 
should  be  met,  which  are  as  follows: 
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(a)  The  displacement  function  must  be  continuous 
within  the  element  and  compatible  between  adjacent  elements. 

(b)  The  function  includes  rigid  body  displacement 
of  the  element. 

(c)  The  function  reduces  an  element  to  a  state  of 
uniform  strain  when  the  element  approaches  a  very  small  size. 

Once  the  shape  and  displacement  function  of  an  ele- 
ment is  chosen,  the  next  important  decision  is  the  choice  of 
material  properties.   The  most  attractive  feature  of  the  fi- 
nite element  method  is  the  versatile  way  of  considering  the 
material  properties.   Theoretically  each  element  can  be 
assigned  different  properties,  which  can  be  changed  in  each 
step  of  analysis.   This  possibility  enables  one  to  analyze 
complex  cases,  like  different  material  properties  at  dif- 
ferent portions  of  a  structure,  nonlinear,  and  anisotropic 
materials,  and  limiting  values  for  certain  element  proper- 
ties, which  permits  accounting  for  slip  or  limiting  the 
magnitude  of  tensile  resistance  between  two  elements. 

The  last  and  final  step  in  the  sequence  of  analysis 
deals  with  techniques  for  numerical  analysis.   Basic  require- 
ments for  success  and  efficiency  of  a  finite  element  analysis 
depends  heavily  on  the  solution  scheme.   Solution  of  a  real 
problem  involves  a  vast  number  of  mathematical  and  numerical 
operations,  which  can  only  be  dealt  with  by  the  use  of  high 
speed  digital  computers.   Hence,  a  good  knowledge  of  computer 
language,  mechanism  and  peripheral  systems  of  computers, 
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numerical  techniques,  data  storage  and  data  structures  are 
essential.  A  few  general  requirements  for  efficient  solu- 
tion are  stated  below. 

(1)  Total  number  of  numerical  operations  should  be 
minimized. 

(2)  Repetitions  of  calculations  should  be  avoided  or 
minimized. 

(3)  Use  of  core-memory  storage  spaces  should  be  as 
low  as  possible. 

(4)  Input  data  should  be  easy  to  prepare  and  store. 
The  output  should  be  easy  to  interpret. 

(5)  Scheme  for  solution  of  stiffness-influence  ma- 
trix should  be  efficient. 

(6)  Ill-conditions  and  "blow-ups"  should  be 
avoided. 

(7)  Solution  should  be  performed  at  the  minimum 
possible  cost. 

3. 3   Selection  of  Finite  Elements  to  Represent  the  Soil. 

Figure  3.1  shows  a  section  through  a  culvert  buried 
in  soil.   The  components,  which  are  believed  to  affect  the 
system's  behavior  are  labeled.   These  are  (1)   the  pipe, 
(2)   soil-structure  boundary  materials,  (3)   bedding, 
(4)   surrounding  zones  of  select  materials  having  prescribed 
states  of  compaction,  (5)   the  earth  fill,  and  (6)  externally 
applied  loads.   To  analyze  a  real  problem  all  of  the  above 
mentioned  components  require  appropriate  consideration. 
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FIGURE    3.1    TYPICAL    SECTION    OF    PIPE    CULVERT. 
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In  addition,  culvert  problems  are  three-dimensional  in 
nature;  in  many  cases,  differentail  displacements  in  the 
longtitudinal  direction  have  adverse  effects  on  performance. 
However,  in  this  study,  emphasis  is  placed  on  better  repre- 
sentation of  soil  materials  and  on  the  soil-pipe  interaction. 
Accordingly,  a  two-dimensional  finite  element  analysis  was 
chosen  to  represent  the  structure. 

Triangular  shaped  elements  were  chosen  to  represent 
the  soil.   As  culverts  usually  have  a  curvilinear  shape 
(often  circular)  triangular  elements  with  one  curved  boundary 
were  selected  (Figure  3.2).   The  element  has  a  total  of  six 
nodes,  three  corner  nodes  and  three  intermediate  nodes  at 
the  midpoints  of  each  side.   The  reason  for  taking  three 
intermediate  nodes  is  given  in  next  section.   To  avoid  ten- 
dencies towards  instability  of  solution,  any  one  side  of  the 
triangle  should  not  be  too  large  compared  to  other  sides, 
i.e.,  the  element  should  not  be  a  'flat'  triangle. 

3. 4   Selection  of  Displacement  Function. 

In  the  finite  element  method  it  is  not  necessary  that 
the  geometry  and  displacements  of  an  element  be  expressed 
by  the  same  ordered  model.   Accordingly  there  can  be  three 
distinct  combinations  between  the  order  of  displacement 
function  and  order  of  geometry.   (1)   Subparametic,  where 
geometry  is  determined  by  a  lower  order  model  than  displace- 
ment, (2)   Superparametric,  in  which  the  reverse  is  true, 
and  (3)   Isoparametric,  where  the  geometry  and  displacements 
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FIGURE  3.2    TRIANGULAR    ISOPARAMETRIC    ELEMENT. 
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are  defined  by  the  same  order  model.   Zienkiewcz,  et.  al. 
(1969) ,  summarized  the  Isoparametric  finite  element  concept, 
a  powerful  generalized  tool  for  formulating  complete  and  con- 
forming elements  with  models  of  any  polynominal  order.   Also, 
Felippa  (196  6) ,  presented  a  rigorous  treatment  of  the  system- 
atic development  of  refined  two-dimensional  elements  for 
plane  stress,  plane  strain,  and  plate  bending.   The  advan- 
tages of  isoparametric  elements  will  be  discussed  later. 

Before  selecting  the  order  of  the  displacement  func- 
tion, one  should  investigate  the  type  of  results  required 
from  the  analysis.   In  the  culvert  problem,  stresses  in  the 
immediate  vicinity  of  the  pipe  are  very  important  factors 
controlling  performance.   Also,  soil  properties  depend 
strongly  on  the  state  of  stress  so  that  regions  of  high 
stress  gradients  must  be  carefully  delineated.   Hence, 
accurate  determination  of  stresses  any  where  within  an  ele- 
ment (including  its  boundaries)  is  required  in  this  study. 
Accordingly,  a  quadratic  displacement  function  was  selected, 
which  ensures  a  linear  variation  of  strain  (and  hence 
stress)  across  the  element. 

Before  treating  the  derivation  of  element  stiffness 
coefficients,  three  additional  factors  had  to  be  discussed. 

(1)  The  natural  coordinate  system 

(2)  Primary  and  secondary  nodes  and 

(3)  Interpolation  displacement  function. 
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The  natural  coordinate  system,  or  "area  coordinates", 
is  a  local  system  defined  for  a  particular  element  which  per- 
mits the  specification  of  a  point  within  the  element  by  a 
set  of  dimensionless  numbers  whose  magnitude  never  exceeds 
unity.   This  system  is  general  and  easy  to  use,  and  integra- 
tion of  element  stiffness  coefficients  is  greatly  simplified. 
In  Figure  3.2,D(x,y)  is  any  floating  point  within  the  element 
(or  at  its  boundary) .   In  the  natural  coordinate  system 
D(£,,  K9r    Kr>)  i    £■!»  £2'  ^3  are  natural  coordinates,  which  are 


defined  as 


An        A2        A3 


51  "  A*  '"  52  "  A*    '    C3  "  A*  (3,1) 

in  which  A*  is  total  area  of  the  curvilinear  triangle  PQR, 
and  A,  is  the  area  of  sub-triangle  DQR,  and  so  on.   The  rela- 
tion between  the  Cartesian  coordinate  system  and  the  natural 
coordinate  system  can  be  determined  as  follows:   (for  conven- 
ience nodal  points  are  designated  by  numbers  given  below  ) 

P  =  1,  Q=2,  R  -  3,  PQ=4,  QR  =  5,  RP=6. 

X  =  q  X;L  +  K2   x2  +  K3   x3 

y  =  h  y±  +  h  y2  +  h  *3  (3-2) 

l  -  2n  =  q    +  C2    +  C3 

in  which 

A, 


and 


2n  =  ^ 


A*  =  A-j^  +  A2  +  A3  +  A4 
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Solving  (3.2),  E,^,  £  ,  £  can  be  determined  in  terms 
of  (x^,  Y^) /  i  =  1,  2,  3  and  D(x,  y)  ,  which  can  be  expressed 
in  a  general  form 

KL   =    (ai  +  b±x  +  CjLy)/2A*,  i  =  1,  2,  3,  ...  (3.2-a) 

where  ai ,  bi  and  c.  are  constants.   In  terms  of  area 
coordinates 


Node  1  =  (1  -  2n,   0,  0) 

Node  2  =  (0,  1  -  2n,  0) 

Node  3  =  (0,   0,   1  -  2n) 

Node  4  =  (|  -  n,  \   -n  ,  0) 

Node  5  =  (  0,  \   ,    j  -    2n) 

Node   6   =   (|,   0,  |  -  2n) 

Primary  (external)  nodes  are  corner  nodes  like  P,  Q, 
R.   Secondary  (internal)  nodes  are  intermediate  midpoint 
nodes  like  PQ,  QR  and  RP. 

By  definition  an  interpolation  function  is  a  function 
which  has  unit  value  at  one  nodal  point  and  zero  values  at 
all  other  nodes.   The  interpolation  functions  are  polynom- 
ials.  The  order  of  polynomial  is  selected  to  satisfy  the 
requirement  of  the  displacement  function.   Let  F.  define  the 
interpolation  displacement  function  for  the  primary  node  1. 
According  to  the  requirement  at  node  1,  F,  =  1  and  F,  =  0 
at  all  other  nodes. 

Let  F,  =  £, (A£,  +  B)  be  a  quadratic  function;  where 
A  and  B  are  constants  to  be  evaluated. 
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At  node  1,  ^  -  1  -  2n ,  F,  =  1 

F1   at  1  -  C1(A(1  -  2n)  +  B)  =  1  (3.3) 

At  node  4  ^  =  X   ~    2n  ,   F,  =  0 

F1  at  4  =  ?1(A  (12~  2n)  +  B)  =   0  (3.4) 

Solving  for  A  and  B  from  (3.3)  and  (3.4) 

A  =  1  -  2n 
B  =   -  1 


So, 

Fi  ■  ^i(r^-2n"^i  -  1}  <3-5) 

For  a  secondary  node  say,  node  4,  F.  =  1  at  node  4 
and  F,  =  0  at  all  other  nodes.   By  inspection,  it  is  seen 
L  =  0  at  nodes  2,  3,  5  and  L  =  0  at  nodes  1,  3,  6. 

Hence,  ^n'^o  =  °  at  nodes  *•»  2'  3'  5'  6 

Let  F.  =  BE,,E,~,   where  B  is  a  constant  to  be  determined. 

At  node  4 ,  g.  =  ^  -r\ ,  C2  =  J  ~   n '   F4  =  1- 

.".  f4  =  B(|-  n)  (|-  n)  -  l 

Solving,   B  =  5- 

(1  -  2n)Z 

So, 

F*  =  o   ci  ro  (3.6) 

4   (l  -  2n)2  51  C2 
Following  similar  procedure  for  all  other  nodes,  the 
interpolation  functions  can  be  developed  and  are  summarized 
below: 
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Fi  =  ^i<T-^-277  ^1  "  1> 

F2  =  ?2(  1  -  4n  ^2  "  1  -  4n) 
F3  =  ?3(2?3  -  1  +  4n) 

F4  =  * J   h  ^2  (3'7) 

4   (l  -  2n)    x 

F5  =  (1  -  4n)    ?2  C3 
*6    (1  -  4n)  ^3  ^1 


In  matrix  notation, 

F  =  iF1      F2   F3   F4   F5   F6)  (3.8) 

F  defines  the  displacement  functions  which  is 
quadratic  in  nature,  normalized  with  respect  to  the  area 
ratio  and  dependent  only  on  the  global  coordinate  system. 

3. 5   Derivation  of  Element  Stiffness  Matrix  for  Linear 
Strain  Triangular  Finite  Element 

Given  the  displacement  functions  and  nodal  displace- 
ments in  both  x  and  y  directions  at  all  nodes,  the  general 
expression  for  displacement  of  point  D  in  Figure  3.2  can  be 

written  as  follows: 

6 

u-E     u.F. 

i=l 

(3.9) 

6 
v-Z    v.F. 

i=l 
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where  u  and  v  are  displacements  in  x  and  y  directions  re- 
spectively and  F.  are  defined  in  (3.7). 

In  matrix  form  (3.9)  can  be  written: 


f  = 


u 


v 


F   0 
0   F 


u 


v 


(3.10) 


in  which  F  as  defined  in  (3.8)  and 


u  =  {u-l   u2   u3   u4   u5   u6} 


v  =  {Vl   v2   v3   v4   v5   v6} 

are  the  nodal  displacement  vectors. 

To  determine  element  strains  the  partial  derivative 
of  equation  (3.10)  is  taken  with  respect  to  x  and  y.   If  the 
vector  e  defines  strain  components,  then  it  can  be  written 
as : 

f    e 


e  = 


r    = 


xyj 


6u 
6x 

i 

r  f  u 
~x 

6v 

6y 

= 

f  V 

~y 

6u 
6y 

+ 

Sv 

6x  * 

f  u 

k  ~y 

•  (3.11) 


~x 


in  which 


6F  _  6F 
~x  '  6x    65- 


6F    6F 
~y  "  6y  "'  65 


65 


6x 

6y 


1    6F 


«5. 


6F 
55- 


6£2     6F 
6x     6? 


65 


fiy 


'3 


2     6F 


65- 


65 


_2 

6x 


65- 


6y 


1    ""     "»2      UJ      "s3 
Taking  derivatives  of  equations  (3.2)   and  (3.7)  and 

multiplying  appropriate  coefficients,  fx  and  f   can  be 
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determined  and  the  expression  for  strains  can  be  written  as 


e  = 


2A* 


Bl 

0 

0 

B2 

B2 

Bl_ 

~4>      o" 

0    ijj 


in  which 


and 


Bn  =  ( 


B.  =  ( 


6B>1  6K2  6£; 

6x  6x  6x 

<S^1  6?2  55; 

6y  6y  6y 


*  = 


ffip  ) 

6-q 

6F 

«52| 

6F 

6C3 

u 


w 


(3.12) 


Equation  (3.12)  can  be  written  in  matrix  form  as: 


e  =  [b] 


where  [b]  = 


2A* 


"Bl   ° 
0    B. 


B2   Bl 


(3.13) 


ij;    0" 
0    ty 


According  to  the  definition  of  stiffness  matrix,  the 
element  stiffness  matrix  k   can  be  written 
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k  = 


[b] 


[D]  ■  [b]  •  d(vol) 


(3.14) 


vol 


in  which  [D]  is  the  elasticity  matrix,  generated  from 
material  properties  and  type  of  analysis  such  as  plane  strain 
or  plane  stress.   Detailed  discussion  about  material  proper- 
ties and  development  of  [D]  matrix  will  be  presented  in  a 
later  section. 

Equation  (3.14)  can  be  expanded  as  follows: 


k  = 


4A' 


• 

' 

J 

t 

^ 

TO. 

T   0 


0   4* 


T  r 

Bi       °        V 


T      T 
B2     B± 


(1) 


(2) 


Dll  D12  D13 


D21  D22  D23 


D31  D32  D33 


Bl 

0 

$ 

<T 

0 

B2 

0 

* 

B2 

Bl 

.  d(vol) 


(3) 


(4) 


(5! 


Multiplying  (2) ,  (3)  and  (4)  yields 


(3.14a) 


k  = 


4A* 


■ 

" 

• 

0 

i 

• 

^cn 

C12 

• 

J 

to: 

0 

^C21 

c 

22__ 

\p         0 
0    ty 


•  d(vol) 


Carrying  out  the  remaining  multiplications  and  for 
constant  thickness  of  element  ' t' 
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k  = 


4A*' 


area 


$   '  Cll  *  t        ^      '    C12  *  *" 


T 

\p     •  c 


21 


*  i>' 


'22 


* 


d(area)  (3.15) 


The  final  step  in  the  evaluation  of  the  stiffness 
matrix  is  to  perform  the  integration  of  all  coefficients  in 
the  matrix  over  the  whole  area  of  the  element.   Each  element 
of  the  stiffness  matrix  is  defined  in  terms  of  area  coordi- 
nates.  To  facilitate  integration  of  these  coefficients  over 
the  area  of  element,  the  following  expression  is  helpful. 


ra  b  _c   ,   ,        a!  b!  c!       _  ,     , 

h    h    h      dx  d*  =  (a  +  b  +  c  +  2)!   2  (area) 


(3.16) 


Area 


Details  of  derivation  of  element  stiffness  matrix 
for  isoparametric  elements,  are  given  in  Appendix  A. 


3. 6   Finite  Element  Representation  for  Pipe  Culverts 

Thorough  review  of  existing  literature  revealed  that 
in  most  cases  finite  elements  used  to  represent  pipes  are 
either  triangular  or  rectangular  in  shape,  and  these  shapes 
were  initially  adopted  for  this  study.   In  the  course  of 
internal  checks  of  the  general  program  the  validity  of  this 
procedure  was  brought  into  question  in  cases  where  the  pipe 
thickness  becomes  small.   The  triangular  element  (with  two 
degrees  of  freedom  at  each  node)  is  incapable  of  considering 
rotation  at  any  node,  which  means  that  flexure  cannot  be 
adequately  represented  if  the  triangle  has  a  linear  strain 
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distribution.   To  examine  this  problem  further,  two  types  of 
problems  were  solved.   In  the  first  type  a  ring  of  uniform 
thickness  having  internal  radius   a,  and  outer  radius   b,  is 
acted  upon  by  two   line  loads  180  degrees  apart.   In  the 
second  problem  relatively  thick  rings  are  acted  upon  by  an 
externally  applied  uniform,  all-around  pressure.   Table  3.1 
shows  the  error  in  the  resulting  change  in  diameter,  Ax% , 
and  of  the  change  in  outer  radius,   Ab% ,  determined  using 
isoparametric   triangular  finite  elements  compared  with  the 
results  of  a  formal  solution  from  the  theory  of  elasticity. 

Table  3.1   Comparison  of  Ring  Solutions 


Line  Loads 

Uniform  Pressure 

a 

b 

b-a 

b-a 
a 

Q. 

O 

Error 

Ax 

a 

b 

b-a 

b-a 

Error 

Ab 

a 

5" 
10" 
10.5" 

6" 
11" 

11.5" 

1" 
1" 
1" 

0.2 
0.1 
0.095 

2.2 
40.0 
60.0 

4" 

10" 

5" 

6" 

14" 

6" 

2" 
4" 
1" 

0.5 
0.4 
0.2 

2.5 

5.0 

11.0 

From  Table  3.1  it  is  seen  that  as  the  ratio  of  thick- 
ness over  internal  radius  decreases,  the  error  increases,  with 
the  effect  being  far  more  pronounced  for  the  line  load  than  for 
the  uniform  pressure ,  which  clearly  indicated  the  need  for  fur- 
ther investigation.   For  this  purpose  a  ring   subjected  to 
two  line  loads  (P)  applied  diametrically  opposite  to  each 
other,  with  internal  radius  a  =  10  inches  and  outer  radii 
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b  =  11/  12,  14,  16,  18,  and  20  inches  was  chosen.   Solutions 
were  obtained  in  three  different  ways: 

(1)  by  elastic  theory 

(2)  using  isoparametric  triangular  finite  elements 
having  linear  strain  distribution 

(3)  using  segments  of  curved  beams  as  finite  ele- 
ments.  The  distribution  of  circumferential  normal  stress 

(a  ).  across  the  pipe  thickness  at  a  section  90  degrees  from 
the  point  of  loading  (spring  line)  is  considered  for  com- 
parison.  In  Figures  3.3  to  3 . 8  the  normalized  circumferential 

°8/2P  r-a 

stress    /  -7-  is  plotted  against  r— —  ratio  (where   r   is  the 

radial  distance  from  center  of  ring) .   The  solid  line  repre- 
sents the  solution  given  by  the  theory  of  elasticity.   Plotted 
points  are  from  finite  element  solutions  using  isoparametric 
triangular  elements,  and  curved  bar  elements.   For  the  latter 
elements  the  stress  distribution  due  to  nodal  moments  is  cal- 
culated including  a  correction  for  element  curvature  (Roark, 
1943) .   From  these  figures  it  is  clearly  seen  that  the  tri- 
angular and  curved  bar  finite  element  solutions  agree  closely 
with  each  other,  and  with  the  formal  solution,  for   

3. 

equal  to  one,     but  for  (b-a/a)  =  0.1  the  use  of  triangular 

elements  introduces  large  errors.   This  is  illustrated  more 

specifically  in  Figure  3.9,  where  the  error  in  afl     is 

max 

plotted  against  the  ratio  of  wall  thickness  to  internal 
radius.   It  is  seen  that  the  error  may  be  unacceptable  when 
the  ratio  of  wall  thickness  to  internal  radius  is  less  than 
about  0 .U .   it  is  concluded  that: 


40 


4  r- 


Thick  Ring  Solution  from  Theory  of  Elosticity 
Triangular    Isoparametric     Element 
Curved- bar    Element 


/  v    / 
\~7 


^^' 


N  =  No.  of  Layers 
of  Triangular 
Elements 


E  =3x10  psi. 
v  =1/3 


0.8 


Normalised     Radial     Distance,  r-^-9 

'  b-a 


FIGURE  3.3     DISTRIBUTION    OF  CIRCUMFERENTIAL 

NORMAL    STRESS   AT  THE    SPRING     LINE.(N=4) 
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Triangular 
Curved-  bar 


-8  - 


0.0 


b  =18  in. 
a  =  10  in. 

b-o 
a 


>0.8 


1 


1 


0.4 
Normalised     Radial     Distance, 


08 

r-  a 
b-a 


FIGURE  3.4    DISTRIBUTION    OF  CIRCUMFERENTIAL 

NORMAL   STRESS  AT  THE    SPRING    LINE.(N-4) 
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1/1 


o 

E 
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|esi|fe  .4 

f  *P 


•  Triangular 
■  Curved-bar 


3 
U 

U 


T3 

a> 

N 

o 

E 
»_ 
o 
z 


■8 


-16 


b=l6  in 
a  =10  in. 
b-a 


0.6 


0.4 
Normalised     Radial     Distance , 


0.8 


r  -  a 
b-a 


FIGURE   3.5    DISTRIBUTION    OF   CIRCUMFERENTIAL 

NORMAL    STRESS    AT    THE    SPRING    LINE.(N-2) 
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FIGURE  3.6   DISTRIBUTION    OF   CIRCUMFERENTIAL 

NORMAL    STRESS  AT  THE    SPRING  LINE.(N-2) 
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60  r- 


—  Thin  Ring  Solution  from  Elastic  Theory, 

Corrected  for  Curvature 
•    Triangular 
■    Curved-bar 


b  =12  in.     b^J  e02 
a  =10  in.      a 


0.4 
Normalised     Radial    Distance , 


OB 


b-a 


FIGURE  3.7    DISTRIBUTION    OF   CIRCUMFERENTIAL 

NORMAL   STRESS   AT  THE    SPRING  LINE. (N- 1) 
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Thin   Ring  Solution  from  Elastic  Theory, 
Corrected  for  Curvature. 

Triangular 

Curved -bar 


0.4 
Normalised    Radial     Distance, 


r-  a 
b-a 


FIGURE    3.8 


DISTRIBUTION    OF    CIRCUMFERENTIAL 
NORMAL    STRESS  AT  THE  SPRING   LINE.  (N- 1) 
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•  N  =  4 
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Ratio    = 


0.4  0.8  1.0 
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Internal  radius 


FIGURE   3.9     ERROR    IN  <remQX    DETERMINED    USING 

ISOPARAMETRIC    TRIANGULAR    ELEMENTS 
VS.  RATIO    OF  PIPE   WALL    THICKNESS 
TO    INTERNAL     RADIUS. 
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(1)  isoparametric  triangular  finite  elements  can  be 

used  to  represent  "thick"  rings,  ( greater  than  o.U)  where 

compression  plays  the  major  role  and  bending  is  of  much 

lesser  importance,  provided  the  ring  thickness  is  subdivided  sufficiently. 

(2)  for  - —  less  than  0.2  bending  becomes  signifi- 
cant and  finite  elements  without  a  degree  of  freedom  in  rota- 
tion can  result  in  significant  errors  in  the  calculated 
stresses,  unless  an  extremely  fine  mesh  is  used. 

(3)  for  ratio  less  than  0.1,  bending  of  the 

3. 

ring  is  dominant,  and  rotation  at  nodes  must  be  considered 
in  the  finite  element  idealization. 

V"»  —  a 

In  the  case  of  corrugated  metal  culverts,  - — 
values  are  nearly  always  less  than  0.1.   In  this  range,  any 
finite  element  using  equivalent  thickness  and  without  rota- 
tional stiffness  at  the  nodes  is  bound  to  have  built-in 
errors  and  triangular  isoparametric  finite  elements  are  not 
suitable  for  use.   This  led  to  the  search  for  a  pipe  element 
which 

(1)  provides  rotational  stiffness  at  the  nodes 

(2)  has  circumferential  and  shear  stiffness 

(3)  can  represent  actual  pipe  thicknesses  that  are 
small  compared  to  the  radius 

(4)  would  not  require  adjacent  soil  elements  to 
have  rotational  stiffness. 

These  requirements  have  been  satisfied  in  a  curved- 
bar    element  which  is  described  in  the  following  section. 
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3 . 7   Curved-Ring  Element 

Martin  (1966)  presented  a  development  for  the  stiff- 
ness matrix  of  curved  beam  elements.   In  this  procedure  a 
section  of  a  circular  ring  is  acted  upon  by  a  normal  force, 
a  tangential  force,  and  by  a  moment  at  each  node.   Each  node 
has  three  degrees  of  freedom  namely,  normal  and  tangential 
displacements,  and  rotation.   Total  moment  (M)  at  any  point 
due  to  applied  forces  and  moments  is  calculated.   Total 
strain  energy   U   is  given  as 


U  = 


2EI 


M2-ds  (3.17) 


where   ds   is  increment  of  length  of  arc  of  bar,  E  is 

the  Young's  modulus  of  pipe  material  and  I  is  the  moment  of 

inertia  of  the  pipe  section. 

Displacements  and  rotations  in  each  degree  of  free- 
dom is  determined  using  Castigliano ' s  theorem  which  yields: 

u  =  Trr-  =   shear  displacement  (3.18) 

oP 

r  tt 

v  =  §77  =   normal  displacement  (3.19) 

oQ 

9  „■  |g  =   rotation  (3.20) 

SM 

where   u   and   P   are  shear  displacement  and  force,  simi- 
larly v   and  Q   are  in  the  normal  direction,  and   9   is 
the  rotation  and  M   is  the  moment.   Substituting  in  the 
equations  for  u,  v  and  6  the  values  of  M,  differentiating 
with  respect  to  P ,  Q  and  M  and  integrating  over  the  arc 
length,  one  obtains  for  node  1 
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u. 


V 


■  [A11]  • 


M 


1/R 


(3.21) 


where  R  is  the  mean  radius  of  arc,  and  A,,  is  a  (3x3) 
matrix  consisting  of  functions  of  known  quantities  like  R, 
EI,  and  3,  the  angle  subtended  by  the  arc  at  the  center. 
Inverting  the  above  matrix-equation 


f  y. 


M, 


f  \ 

Ul 


-  [Kn]  '  * 


V. 


(3.22) 


where  [K,-^  =  inverse  of  [A.^] 

For  the  other  and  end  of  the  arc,  (node  2) ,  another 
matrix  equation  can  be  written  in  a  similar  fashion. 


\Qn\      =  [K22l 


"2 


u. 


(3.23) 


Considering  total  equilibrium  of  the  arc  under  the 
system  of  forces  and  moments, 
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'    1 

P2 

Ul 

Q2 

\       =     [K21]   •   - 

Vl 

M2 

61 

(3.24) 


Combining  all  these  matrix  equations 


ul 

vl 

■  01 

U2 

V2 

02 

^11         k12 


k21      k22 


6x6 


6x1 


'V 

Ql 

Ml 

P2 

Q2 

M2 

(3.25) 


6x1 


where  k, _  =  transpose  of  matrix  k„, .   The  result  is  a 
(6x6)  stiffness  matrix  for  a  curved  beam  element.   Details 
of  this  derivation  can  be  found  in  Martin,  (1966)  and  are 
given  in  Appendix  B.   The  final  step  in  the  formulation  of 
the  stiffness  matrix  for  a  curved  beam  element  is  the  trans- 
formation to  the  global  coordinate  system,  which  can  be 
accomplished  by  matrix  multiplication  of  the  stiffness  matrix 
with  the  transformation  matrix  whose  coefficients  consist  of 
appropriate  direction  cosine   (For  details  see  Appendix  B) . 
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3. 8   Soil-Pipe  Interaction  Element 

Many  researchers  have  expressed  the  need  for  consid- 
ering interaction  between  pipe  and  enclosing  soil  medium. 
This  interaction  can  involve  slip  between  the  pipe  and  the 
soil  at  shear  stresses  less  than  the  shear  strength  of  the 
soil,  especially  if  the  pipe  is  wrapped  in  a  polysterene 
film.   A  suitable  interaction  element  should  have  the  follow- 
ing features: 

(1)  In  contact  with  granular  soils,  no  net  tensile 
resistance  in  the  normal  direction 

(2)  Negligible  displacement  in  the  normal  direction 

(3)  Small  shear  displacements  as  long  as  the  ratio 
of  shear  stress  to  normal  stress  is  less  than  a  given  lim- 
iting value  (limiting  wall  friction) 

(4)  When  the  ratio  of  shear  stress  to  normal  stress 
exceeds  the  limiting  value,  the  shear  stress  is  sustained  but 
there  is  negligible  resistance  to  shear  displacement. 

Accordingly,  an  element  of  zero  thickness  was  chosen, 
which  is  capable  of  offering  resistance  to  shear  deformation 
until  a  limiting  ratio  of  shear  to  normal  stress  is  exceeded. 

Goodman,  Taylor  and  Brekke  (19  68)  presented  a  model 
for  the  mechanics  of  jointed  rock.   Discontinuities  or 
joints  in  the  rock  mass  are  considerd  as  different  material 
with  distinct  properties  which  are  different  from  those  of 
adjoining  rock  blocks.   To  represent  joints  in  terms  of 
finite  elements,  the  authors  derived  the  stiffness  matrix  for 
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these  joints  based  on  the  minimum  potential  energy  concept. 
The  procedure  presented  by  Goodman  et.  al.  (1968)  had  been 
used  in  this  study  to  develop  the  stiffness  matrix  for 
interaction  elements.  The  element  has  four  nodes,  each  node 
having  two  degrees  of  freedom  in  the  normal  and  tangential 
direction  of  the  interaction  element.   A  pair  of  nodes  of  the 
interaction  element  (Figure  C-l)  separated  by  a  finite  dis- 
tance are  connected  to  two  consecutive  nodes  of  a  pipe  seg- 
ment.  The  other  pair  of  nodes  are  connected  to  the  adja- 
cent soil  elements.   Initially  the  normal  distance  between 
these  two  pairs  of  nodes  is  zero.   Properties  for  the  inter- 
action element  are  described  in  Chapter  IV,  and  detailed 
derivation  of  the  stiffness  matrix  for  this  element  with 
eight  degrees  of  freedom  is  presented  in  Appendix  C. 
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CHAPTER  IV 
MATERIAL  PROPERTIES 

4 . 1   Introduction 

One  of  the  primary  advantages  of  the  finite  element 
method  is  the  capacity  of  handling  realistic  material  proper- 
ties which  may  differ  from  place  to  place  in  the  structure. 
Referring  to  Figure  3.1,  there  are  three  basic  types  of 
material  behavior  in  the  culvert  problem;  (1)   pipe, 
(2)   soil-pipe  interaction  and  (3)  soil  medium.   Depending 
upon  the  number  of  different  compacted  zones  several  soil 
types  with  distinct  properties  may  exist.   At  the  interface 
between  the  pipe  and  soil  the  behavior  is  different  from 
the  properties  of  either  of  the  two  media.   This  will  be 
modeled  by  an  interaction  element  of  zero  thickness  with 
its  own  inherent  properties. 

It  is  well-known  that  an  isotropic  elastic  material 
requires  two  independent  parameters  to  define  its  proper- 
ties.  Most  commonly  used  elastic  parameters  are  Young's 
modulus  (E)  and  Poisson's  ratio  (v) .   For  linearly  elastic 
material  both  E  and  v  remain  constant  over  the  entire  range 
of  loading  and  unloading.   Soil  is  not  an  ideal  material  in 
this  respect;  on  the  contrary  its  stress-strain  relation  is 
non-linear  and  stress  path  dependent.   Natural  clay  soils, 
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up  to  a  certain  limit,  show  linearly  elastic  behavior,  but 
sands  are  nonlinear  and  inelastic  even  at  low  stress  levels. 
For  this  reason,  analyses  that  assume  the  soil  to  be  a 
linearly-elastic  material  are  not  appropriate. 

4 . 2  Pipe  Properties 

Commonly  used  culvert  materials,  such  as  steel, 
aluminum  or  concrete  show  linear  elastic  behavior  both  in 
compression  and  in  tension   up   to   the   yield   point. 
Within  this  stress  range  it  is  reasonable  to  assume  that  the 
properties  of  the  pipe  materials  remain  constant.   It  is 
recognized  that  before  collapse  concrete  pipes  will  crack 
and  steel  culverts  will  yield.   However,  in  this  study  it 
was  decided  to  concentrate  on  the  behavior  of  the  soil 
medium  and  the  soil-pipe  interaction  effects.   Accordingly, 
it  was  initially  assumed  that  values  of  E  and  v  for  the  pipe 
would  remain  constant  and  to  consider  only  problems  in  which 
the  yield  point  stress  of  the  pipe  materials  is  not  exceeded. 

4. 3  Soil-Pipe  Interaction 

Introduction 

A  number  of  investigators  conducted  model  and  field 
tests  on  flexible  and  rigid  pipes  buried  in  soil;  Krizek,  et, 
al.  (1974);  Kirkland  and  Walker,  (1972);  Neilson,  (1972); 
Brown,  (1966);  Davis,  (1966,  1969);  Abel,  (1973)  and 
Allgood,  (1972).   Most  of  this  research  shows  that  measured 
results  deviate  from  the  analytical  predictions  even  under 
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controlled  test  conditions.   Among  the  factors  considered  to 
explain  the  discrepancy  is  the  interaction  between  soil  and 
pipe,  Krizek  (1972)  ,  Linger  (1972) ,  Corotis  et.  al.  (1974) , 
Abel  et.  al.  (1973)  and  others.   The  term  soil-structure 
interaction  refers  to  the  general  phenomenon  involved  in  the 
behavior  of  a  buried  structure  as  a  result  of  different  prop- 
erties at  the  interface  of  the  structure  and  the  surrounding 
soil  medium  in  response  to  the  imposed  loading  system.   This 
concept  has  existed  in  the  literature  for  a  long  time,  how- 
ever, until  recently,  the  phenomenon  received  little  atten- 
tion.  Linger  (1972)  presented  a  historical  development  of 
soil-structure  interaction.   The  concept  may  be  stated  as 
the  process  of  pressure  distribution  between  soil  and  struc- 
ture, which  depends  on  the  degree  to  which  relative  dis- 
placements along  the  interface  has  mobilized  the  contact 
shear  strength. 

In  the  conduit  problem,  the  existence  of  soil- 
structure  interaction  was  first  recognized  by  Spangler  and 
Mars ton  (1941)  and  they  introduced  the  concept  of  the  "im- 
perfect ditch  method."   As  a  result,  different  methods  of 
pipe  installation  have  been  developed.   For  example,  in  the 
"Soil-Arch"  method  a  layer  of  material  is  compacted  around 
the  pipe  after  a  relatively  stiffer  material  has  been 
placed  at  both  abutments.   The  soil  arch  helps  transfer 
load  from  overburden  to  a  base   away   from  the  pipe. 
Fisher  (1969)  described  the  installation  of  "super-span" 
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culverts  whose  successful  performance  he  attributes  to  slip 
along  the  soil-pipe  interface  and  to  transfer  of  load  by 
arching  to  well-compacted  zones  beside  the  base  of  the  arch. 

Properties  for  the  Interaction  Element 

Interaction  between  soil  and  structure  may  be  looked 
at  as  a  phenomenon  which  occurs  as  a  result  of  existence  of 
one  material  next  to  another.   Because  of  the  fact  that  the 
two  materials  have  different  properties  the  behavior  at  the 
interface  is  different  from  that  of  either  material  by  it- 
self.  The  mechanism  of  shear  transfer  between  a  structure 
and  clay  soil  is  little  understood  whereas  for  sandy  soil 
some  information  is  available  (Brumund  and  Leonards,  (1973)). 
Interface  friction  depends  on  several  factors,  among  which 
are  the  normal  pressure,  roughness  of  contact  surfaces, 
relative  displacements,  and  environmental  conditions.   For 
a  granular  soil,  the  following  characteristics  are  recognized 
as  representative  of  behavior: 

(1)  In  two-dimensional  representation,  the  inter- 
action element  should  approach  zero  thickness. 

(2)  Tension  normal  to  the  pipe  cannot  be  tolerated. 

(3)  Compression  displacements  normal  to  the  pipe 
are  negligible. 

(4)  The  shearing  resistance  is  purely  frictional  in 
nature,  and  is  proportional  to  the  normal  stress. 

(5)  The  resistance  to  shear  deformation  is  high 
until  the  ratio  of  shear  stress  to  normal  stress  exceeds  a 
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limiting  value.   Beyond  this  limit  ,   the   shear  stress  is 

maintained  but  full  slip  is  permitted. 

Based  on  the  above  mentioned  features,  the  model  for 

the  properties  of  the  interaction  element  is  shown  in  Figure 

4.1.   High  values  for  stiffness  in  the  normal  and  shear 

direction  (En  and  E   in  equation  C.9)  are  initially  assigned, 

If  after  application  of  an  increment  of  load,  the  net  normal 

stress  is  tensile,  E   and  E   are  reduced  to  zero.   If  the 

n      s 

net  normal  stress  is  compressive  but  the  stress  ratio  ex- 
ceeds the  limiting  value,  then  E   is  reduced  to  zero  while 
En  is  left  unchanged,  which  simulates  slip  in  interaction 
element. 

4 . 4   Soil  Properties 

In  general,  stress-strain  properties  for  soil  are 
non-linear  in  nature;  even  at  small  strains,  the  stress- 
strain  curves  are  slightly  nonlinear  and  there  are  plastic 
components  of  strain.   However,  the  amount  of  plastic  strain 
is  not  large  in  the  initial  stages.   Frydman  and  Zeitlen 
(1969)  concluded  from  triaxial  tests  performed  on  dense 
dune  sand  and  glass  microspheres  that  only  insignificant 
amounts  of  plastic  strain  occur  before  a  yield  point  is 
reached.   A  major  factor  which  influences  the  stress-strain 
behavior  for  sand  is  the  confining  pressure  (Lee  and  Seed, 
1967,   Duncan  and  Chang ,  (1970)  . 

Three  stress-strain  models  for  soils  are  commonly  in 
use:   linear  elastic,  elastoplastic  and  non-linear.   A 
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Same  coordinates, 
(zero  thickness, 
after  Goodman, 
Taylor  S   Brekke, 
1968) 


Es  =    En    until     r/cr   =  tan  8  (interaction    friction) 
then,  Eg**0 


FIGURE  4.1  PROPERTIES     OF    INTERACTION 

ELEMENT. 
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drawback  of  the  linear  elastic  model  is  that  it  does  not  con- 
sider shear  dilatancy.   The  state  of  stress  at  a  point  can 
be  expressed  in  terms  of  octahedral  normal  stress 

Oct    J    1     2     3 
and  the  octahedral  shear  stress : 


Toct  =  J  /  (al  "  a2)2  +  (a2  "  a3)2  +  (a3  "  °1)2 

(4.2) 
Frydman  and  Zeitlen  (1969)  showed  that  in  triaxial 
tests  with  a      .    constant,  no  volume  change  took  place  before 

a  particular  ratio  of  t  ./a      .    was  exceeded  regardless  of 

*  oct   oct  3 

the  value  of  a   ..   The  point  at  which  volume  change  started 
oct       c 

was  called  the  yield  point,  after  which  the  stress-strain 
curve  showed  increasing  change  in  curvature;  this  aspect  is 
important  for  determination  of  deformations. 

Several  models  for  simulating  nonlinear  stress-strain 
behavior  of  soils  are  available  in  literature.   Chen  (1948) 
concluded  from  triaxial  tests  on  cohesionless  soils  that 
axial  strain  can  be  expressed  in  terms  of  an  exponential 
function  of  deviator  stress  for  a  stress  range  of  15%  to  50% 
of  maximum  deviator  stress.   A  hyperbolic  function  was  pro- 
posed by  Kondner  (1963)  to  represent  stress-strain  relations 
for  cohesive  soils.   He  used 

(a,  -  a,)  =    I   .  (4.3) 

1     3     a  +  be 
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in  which  'a'  and  'b'  are  constants  to  be  evaluated  from  ex- 
perimental results.   Taking  derivatives  of  equation  (4.3) 
one  can  find  the  value  of  tangent  modulus.   He  also  proposed 
methods  to  determine  the  constants  'a'  and  'b'. 

Brinch  Hansen  (1963)  proposed  a  parabolic  relation- 


ship  as  a   =   .   Where   r   and   s   are  constants  of 

c  r  -  s«  e 

the  parabola.   The  initial  tangent  modulus  value  is  indefin- 
ite at  e  =  0 . 

Duncan  and  Chang  (1970)  used  the  hyperbolic  model 
for  stress-strain  relations  proposed  by  Kondner  (1963)  modi- 
fied to  account  for  the  stress-dependency  of  the  model.   They 
also  developed  relations  for  tangent  modulus  which  is  suit- 
able for  use  in  incremental  analysis,  such  as  the  finite 
element  analysis.   The  main  advantage  of  using  the  tangent 
modulus  instead  of  initial  modulus  is  that  in-situ  stresses 
can  be  taken  into  account  and  incremental  analysis  can  be 
performed,  taking  state  of  stress  and  strain  dependency 
effects  into  account.   Their  expression  for  tangent  modulus 
was  given  as 

n 


Et   = 


Rf    (1   -   Sin    <J>)     (a1   -   a3) 
2   C      •    Cos    <$>   +      2a3    •    Sin   <jj 


a3 
K  •  p   •  (— ) 
^a    p 
^a 


(4.4) 


in  which  R_  =  failure  ratio,  c  and  <|>  are  cohesion  and  fric- 
tion angle  values  of  the  soil,  p   is  the  atmospheric  pres- 
sure,  K  and  n  are  constants  to  be  determined.   This  model  is 
useful  for  soils  which  show  hyperbolic  stress-strain 
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relations.   For  soils  in  which  the  stress-strain  relation 
deviates  from  hyperbolic  shape,  use  of  the  model  introduces 
significant  errors.   Duncan  and  Chang  (1970)  used  the  pro- 
posed relationship  for  tangent  modulus  in  calculation  of 
footing  behavior,  which  show  significant  improvements  on 
previous  analytical  predictions.   Domaschuk  and  Wade  (1969) 
proposed  the  use  of  bulk  modulus,  K  and  shear  modulus  G,  as 
opposed  to  Young's  modulus,  E  and  Poisson's  ratio,  v  be- 
cause K  and  G  can  be  evaluated  independently. 

Lade  (1972)  did  an  extensive  study  of  stress-strain 
relations  for  soil.   He  presented  a  three-dimensional, 
elastoplastic  stress-strain  relationship  for  cohesionless 
soils.   This  development  is  based  on  the  concept  that  during 
loading  both  elastic  and  plastic  strain  occur.   The  elastic 
part  is  evaluated  using  elastic  theory  and  plastic  component 
is  determined  using  plastic  stress-strain  relations. 

e     p 

Ae   =  Ae   +  Ae^  , .  c , 

xxx  (4.5) 

e 

in  which  Ae   is  the  total  increment  of  strain,  Ae   is  the 
x  x 

elastic  and  Ae   is  the  plastic  component  of  strain.   This 
model  is  based  on  three  requirements: 

(1)   There  must  exist  a  yield  surface  such  that  if 
the  stress  condition  is  within  this  surface,  the  soil  will 
behave  elastically.   If  the  point  is  outside  the  surface,  it 
will  deform  plastically.   When  the  point  is  on  the  surface, 
it  represents  yield  a  condition. 
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(2)  A  flow  rule  should  relate  plastic  stresses  to 
plastic  strain. 

(3)  A  work-hardening  law  is  required  from  which  the 
magnitudes  of  strain  increments  caused  by  a  given  stress 
increment  can  be  determined.   The  yield  criterion  adopted  is 
of  the  form: 

Xl   "  kl  I3  =  °  (4'6) 

f  =  I13/I3 

in  which  I,  and  I3  and  the  first  and  third  principal  stress 
invariants  and  k,  in  a  constant  whose  values  depend  on  the 
type  of  sand,  and  f  is  stress  level.   When  an  increment  of 
stress  is  applied,  the  material  will  yield  plastically  if 
df  >  0,  and  will  yield  elastically  if  df  =  0  or  df   <  0. 
Lade  (1972)  assumed  the  following  function  for  plastic 
potential. 

g  -  Ix3  -  k2  I3  (4.7) 

in  which  g  is  the  plastic  potential  and  k„  is  a  constant  for 
given  value  of  f.   The  plastic  stress-strain  relation  is 
given  by : 

d  ep.  .  =  dX  •  -^2_  (4.8) 

il  do .  . 

in  which  de  . .  is  the  plastic  strain  increment  in  i ,  j 

directions  due  to  an  increase  in  stresses  da . .  in  the  same 

iD 

directions;  dX    is  constant,  which  relates  to  the  work  harden- 
ing rule.   Total  plastic  work  done  due  to  applied  stress 
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system  can  be  written  as  a  function  of  f. 


w  =  F(f) 
P 

wp  =  plastic  work,  f  =  I   /I 


(4.9) 


Lade  (1972)  has  shown  that  for  cubical  stress  conditions 
dW 


AA  = 


3g 


(4.10) 


Taking  derivatives  for  all  strain  components  in 
equation  (4.8)  substituting  for  I„  and  AX  the  plastic  strain 
components  can  be  written  as: 


AeJ 

Ae 
Ae 
Ae 


x 

P 

y 
p 

z 

p 

-yz 

AeP 

Ae 


zx 

P 


xy 


K, 


dW 


3(i: 


K2Il) 


3  T2 


a      +  t 
y   z     yz 


3   2 

=■1,-0   a   +  t" 

K„  1     z   x     zx 


2 

z: 

2 


3   2 

K2  1    x   y    xy 

2a   t    -  2  t    t 
x  yz      xy   zx 


2a   t    -  2  x    x 
y   zx      xy  yz 


2a   x    -  2  x    t 
z   xy      yz   zx 


(4.11) 
Lade  (1972)  presented  methods  for  evaluating  the  constants 
of  equation  4.11  from  standard  tests.   He  used  these  rela- 
tions in  a  variety  of  cases  including  cubical  triaxial  tests 
and  has  shown  that  predicted  results  compare  well  with  the 
experimental  results.   To  investigate  the  influence  of 
intermediate  principal  stress,  a_,  Bishop  (1966)  proposed  a 
parameter  b,  which  is  defined  as 
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°2    -  °3 
b  =  -= (4.12) 

al  -  a3 

Lade  (1972)  presented  results  for  range  in  values  of   b 
varying  from  0  to  1 . 

Need  for  Incremental  Analysis 

Geotechnical  engineering  construction  often  involves 
either  filling  or  making  excavations  in  soil  which  is  done  in 
layers.   Classical  methods  of  analysis  do  not  take  into 
account  this  real  situation.   For  example,  in  classical 
analysis  in  the  case  of  construction  of  an  embankment,  a 
decision  has  to  be  taken  regarding  full  shear,  full  slip  or 
partial  slip  between  layers.   In  such  cases  use  of  finite 
element  method  has  an  important  advantage  over  many  other 
existing  methods. 

Goodman  and  Brown  (1963)  were  among  the  first  to 
recognize  the  importance  of  considering  the  effects  of  con- 
struction or  excavation  in  layers.   If  construction  is  ac- 
complished in  layers  and  the  soil  possesses  nonlinear  be- 
havior, simulation  of  this  process  in  analytical  models  is 
important.   Even  if  the  material  is  linearly  elastic  and 
deformations  are  small,  construction  in  layers  would  have  an 
effect  because  of  incompatibility  at  the  interface  between 
layers.   One  of  the  great  advantages  of  the  finite  element 
method  is  its  ability  to  simulate  the  layer-by-layer  con- 
struction sequence.   To  minimize  the  cost  of  such  analyses, 
the  general  practice  is  to  use  actual  lift  heights  in  the 
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zones  of  major  influence  on  the  structure  and  relatively 
larger  lift  heights  elsewhere. 

The  influence  of  the  construction  sequence  is  of 
importance  in  many  other  types  of  problems;  for  example, 
retaining  walls,  locks,  dams  and  buried  conduits.   In  par- 
ticular, for  problems  of  flexible  culverts  in  soil,  many 
investigators  have  pointed  out  that  slip  between  pipe  and 
soil  may  occur  and  this  might  have  an  important  effect  on 
the  behavior  of  the  system.   To  account  for  this  slippage, 
the  characteristics  of  which  depend  on  normal  thrust  on  the 
interface,  incremental  finite  element  analyses  is  very 
convenient. 

The  construction  of  a  culvert  may  be  considered  as  a 
three  phase  process.   The  first  phase,  a  trench  may  be 
excavated  (changing  the  free-field  stress  system)  and  the 
bedding  is  prepared.   In  the  second  phase,  the  pipe  is 
placed  and  initial  back  packing  takes  place  up  to  the  crown. 
During  this  phase  the  pipe  is  not  fully  constrained  by  sur- 
rounding soil  and  considerable  displacements  may  occur.   The 
third  phase  includes  filling  and  backpacking  operations  up 
to  the  top  of  the  embankment  or  trench.   These  three  pro- 
cesses take  place  step-by-step  and  may  not  be  carried  out 
symmetrically  with  respect  to  the  pipe.   Accordingly,  con- 
struction effects  may  be  simulated  by  incremental  analysis 
but  the  real  situation  can  not  be  fully  recreated. 
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In  this  study  a  combination  of  incremental  and  iter- 
ative procedures  has  been  used  depending  upon  the  require- 
ments of  the  problem.  In  general  the  incremental  approach 
is  used  in  the  analysis  until  some  ill-condition  like  ten- 
sion in  the  soil  or  slip  in  the  interaction  element  occurs. 
In  such  cases  an  iterative  approach  must  be  used.  Details 
of  this  approach  are  presented  in  a  later  chapter. 

4 . 5   Method  of  Considering  Stress-Strain  Parameters 

For  incremental  finite  element  analyses,  the  follow- 
ing factors  are  important  in  selection  of  soil  parameters: 

(1)  Stress-strain  relations  are  non-linear  and 
dependent  on  the  state  of  stress  and  the  stress  path. 

(2)  A  tangent  modulus  and  tangent  Poisson's  ratio 
is  more  representative  of  non-linear  behavior 

(3)  A  single  stress-strain  formulation,  which  is 
reasonably  accurate  over  a  wide  range  of  conditions,  prob- 
ably does  not  exist. 

Solutions  of  nonlinear  problems  using  the  finite 
element  method  can  be  done  in  three  ways  -  (1)   stepwise  or 
incremental  procedure,  (2)   iterative  technique  and  (3)   com- 
bination of  incremental  and  iterative  methods.   In  the 
incremental  method  the  total  load  is  divided  into  predeter- 
mined numbers  of  increments.   Results  such  as  stresses,  dis- 
placements of  each  step  are  added  to  the  results  at  the  end 
of  previous  step.   This  process  is  continued  till  the  total 
load  is  applied.   The  underlying  assumption  in  this  method 
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is  that  the  material  behaves  linearly  during  each  increment. 
The  accuracy  of  this  approach  depends  largely  on  the  number 
of  increments.   The  source  of  error  is  illustrated  in 
Figure  4.2. 

In  Figure  4.2,  A  and  B  are  two  points  representing 
an  increment  on  the  stress-strain  curve.   AC  is  the  tangent 
modulus  at  point  A,  which  is  different  from  the  tangent 
modulus  at  point  B.   Use  of  a  tangent  modulus  value  either 
at  A  or  at  B  for  the  increment  will  introduce  an  unavoidable 
error.   Somewhat  better  result  may  be  expected  by  using  an 
average  value  of  the  tangent  modulus  between  points  A  and  B, 
(EF) .  But  there  are  two  difficulties:   (1)   before  analyzing 
for  the  increment,  point  B  is  undefined  and  (2)   segment  of 
the  curve  between  A  and  B  may  be  strongly  curved.   Using  an 
iterative  technique   the  error  can  be  reduced.   One  drawback 
of  this  approach  is  that  the  solution  becomes  expensive.   To 
eliminate  these  sources  of  errors,  cubic  spline  functions 
were  used  to  simulate  the  actual  stress-strain  curve.   The 
basic  principles  of  spline  functions  are  consistent  with  the 
discretization  principle  of  finite  elements. 

Desai  (1971)  presented  a  method  for  using  spline 
functions  for  nonlinear  finite  element  analysis.   From  a 
given  set  of  stress-strain  curves  for  different  confining 
pressures,  another  set  of  splines  can  be  generated  for  tan- 
gent modulus  versus  confining  pressure  for  any  given  strain 
or  stress  ratio  as  desired.   Desai  (1971)  solved  a  footing 
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Strain 


FIGURE  4.2      ERRORS    IN    INCREMENTAL    ANALYSIS. 
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problem  with  nonlinear  stress-strain  properties  using  spline 
functions,  keeping  Poisson's  ratio  constant,  which  shows 
better  agreement  than  assuming  hyperbolic  stress-strain  rela- 
tions.  Desai  (1972)  presented  another  method  of  approximat- 
ing stress-strain  relationships  using  bicubic  spline  func- 
tions.  In  this  process  continuity  in  all  three  spaces 
(e.g.,  stress,  strains  and  confining  pressure)  is  ensured. 
However,  due  to  computational  difficulties  of  bicubical 
spline  functions,  cubical  spline  functions  were  used  in  this 
study.   The  approach  taken  is  slightly  different  from  that 
proposed  by  Desai  (1971) .   The  advantage  of  using  spline 
functions  in  incremental  analysis  is  that  the  parameter  in 
question  can  vary  continuously  within  a  stress  (or  strain) 
interval. 

Spline  Functions 

In  this  study  a  general  method  for  handling  incre- 
mental loadings  and  non-linear  material  properties  has  been 
adopted  making  use  of  spline  functions.   Ahlberg  et.  al. 
(1967)  presented  a  comprehensive  mathematical  background  of 
a  variety  of  spline  functions  and  their  application.   In 
Figure  4.3,  (a,b)  represents  an  interval  A  of  a  typical 
stress-strain  curve.   The  interval  has  been  subdivided  into 
a  number  of  subintervals  so  that 

A  :  a  =  x   <  xn  <  . . .x   =  b 
o     1       N 

and  corresponding  ordinates  are: 

y  :  y  ,  yv   y2,  .  .  .  yN  (4.13) 
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in  which  Y  is  the  vector  of  ordinate  at  each  x.  value. 

1 

Spline  function  may  be  defined  as 

SA  (x)  =  Y  (4.14) 

S   is  continuous  together  with  its  first  and  second 

derivatives  in  the  interval  (a,  b)  and  coincides  with  a  cubic 

in  each  subinterval,  i.e., 

x .  ,  <  x  <  x .  for  j  =  1 ,  2 ,  . . .  N  and 

(4.15) 

S.(x.)  =  Y.  for  j  =  0,  1,  2,  ...N 
A   j      j      J 


The  second  derivative  of  S.  may  be  considered  as 
"moment"  of  spline  (M . ) . 

x  . -x       x-x._, 
S"  (x)  =  M.  1   -£ —  +  M,    hJ  (4.16) 

3  j  J 

where  h.  =  x.-x.  ,  is  the  size  of  subinterval  and  x  is  a 
D     J      D-l 

variable  point.   Integrating  equation  (4.16)  twice, 
M.  n  (x-x.)3    M.    (x-x.   )3 

sa(x)  =  -  -Jt — r1-  +  h1 e1^-  +  cix  +  c2  (4-17) 

D  D 

in  which  C,  and  C?   are  constants  to  be  determined  from  the 
end  conditions  at  point  j-1  and  j  which  are: 


at    x  =  x.^,    SA(xj_1)  =  Y.^ 

and   x  =  x.   ,    SA  (x  )  =  Y 

After  determining  C,  and  C„  and  substitution  in 
equation  (4.17)  yields: 
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SA(X)  =  6^  (XJ"X)   + 
3     J 


6h.  (x-x.   ) 


M .  ,  h  .         x . -x 


M.h. 
+  (Y.   -    I    ]   )  •  ( 


x-x 


h. 
D 


izi 


) 


(4.18) 


To  ensure  continuity,  the  first  derivative  of 
equation  (4.18)  is  taken  and  the  following  continuity 
condition  is  applied   (Figure  4.4). 


S!  (x.-dx)  =  _  .     S* . (x.+dx) 
Lim    A]       Lim     A   j 


dx->-o 


dx->o 


(4.19) 


After  simplification, 

h,  +  h 
-i- 
j-l 


N  M. 
6 


j+l 


3         ]     6 

h. ,  ,  h .  i 

D+l  1  J 


(4.20) 


Similar  equations  are  written  for  each  interval  and 
can  be  written  in  a  matrix  form  as  below 


2(h1+h2) 


2(h2+h3) 


0 
h. 


2(h3+h4) 


•          •          • 

Ml 

dl 

.    .    . 

M2 

d2 

.    .    . 

• 

M 
2 

= 

d3 

2(h  +h    ... 

n      n+1) 

M 
n 

d 
n 

(4.21) 
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FIGURE  4.4    SUBINTERVAL    OF    SPLINE     FUNCTION. 
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In  equation  (4.21)  the  vector  M.  is  unknown,  but 
can  be  solved  for    using  matrix  algebra.   For  any  value  of 
x  in  x.  i  <  x  <  x.,  the  ordinate  can  be  determined  using 
equation  (4.18).   But  the  primary  interest  is  to  find  the 
slope  at  any  point  which  in  fact  is  the  tangent  modulus. 
Slope  of  a  cubic  spline  is  determined  by  taking  the  first 
derivative  of  equation  (4.18),  which  can  be  written  as: 

M._        0        M.  ~ 

sk  (x)  =  -  21P  (xrx)   +  2h-  (x-xj-i> 

Y.-Y.        M.-M.  , 
+   _J 2—L      -      _J J L   h  (4  22) 

D 

In  this  study,  taking  first  derivative  as  in 
equation  (4.22)  is  avoided  by  starting  with  slopes  rather 
than  ordinate  values. 

Say  S!"  (x.)  =  m.  =  slope  at  point  j.   In  the  inter- 
val (x.  .  ,  x.)  the  expression  for  S  (x)  can  be  written  as: 
j-1     j  m 

2  2 

(x.-x)   (x-x.  ,  )        (x-x...)  (x.-x) 

S(x)  =  m.  ,  —J * 2^-  -  m.  i-± J 

m       ^_1       h2  D       h2 

D  D 

(x-x)2  [2  (x-x    )  +  h.] 

+  Y.  ,   — J _ 2_J: ^_ 

3"1  h3 


(x-x.  ,)2[2(x.-x)  +  h  ] 
+  Y.     3_J: 1 L_  (4.23) 

=>  h3 

D 

and  slope  s' (x)  is  given  by 

c        m 
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(x.-x)  (2  x.'   +  x.-3x) 
S'(x)  =  m.  -   — 3 lli L 


m 


j-l 


(x-x._.) (2  x.  +  x._,  -  3x) 
_  m  .  3 J 2 


+  6  •  1    3->         •  (X.-X)  •  (X-X.  x)    (4.23a) 


Taking  second  derivative  of  equation  (4.23),  en- 
forcing continuity  condition  as  x  approaches  x .  from  both 
sides,  it  can  be  written  in  matrix  form. 


2        \i  0 

o 


A1      2 


0         A. 


2      U. 


ln-l    2      y 


n-1 


A         2 
n 


m 


hk 


m. 


m 


n 


n 


(4.24) 


in  which 


A  .    = 


J±l 


j        h.+h... 


and  »»  .    =   1   -    X . 


and 


C.    =   3A .         J    u  +      3y .    • 
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Whence  vector  m  can  be  solved  for.      Using  appropriate 

values  of  in,  the  slope  at  any  point  can  be  determined  by 

equation  (4.23a). 

In  Figure  4.5(a)  typical  stress-strain  curves  for 

different  confining  pressures  (a-J  are  shown.   Tangent  moduli 

for  nonlinear  material  behavior  can  be  determined  using 

spline  functions  as  follows. 

(1)   From  the  given  data  points, A.,  u.  and  C.  are  determined 

and  the  tridiagonal  matrix  along  with  right  hand  side  of 

equation  (4.24)  is  formed.   The  system  of  equations  are 

solved  for  values  of  m.   (2)   In  the  next  step,  slopes  at 

the  curve  (which  is  tangent  modulus)  at  predetermined  failure 

stress  ratios  (  x   ,/t     )  are  evaluated  using  equation 

oct'  oct  '  3   * 

4.23.   Steps  (1)  and  (2)    are  repeated  for  all  curves. 
(3)   Another  set  of  splines  are  used  to  fit  tangent  modulus 
and  confining  pressure  (a   . )   for  various  values  of  failure 
stress  ratios.   Because  only  values  of  ordinates  are  re- 
quired, equation  (4.21)  is  used  in  this  stage.   Values  of 
tangent  modulus  for  a  given  stress  ratio  and  for  an  interval 
of  confining  pressure  can  be  determined  by  using  equation 
(4.18).   Figure  4.5(b)  shows  curves  for  tangent  modulus 
versus  confining  pressure  for  different  octahedral  shear 
stress  ratios  obtained  from  the  stress-strain  curves  shown 
in  Figure  4.5a.   Using  this  process,  finite  element  analysis 
can  be  performed  for  nonlinear  stress-strain  relations, 
which  in  turn  are  influenced  by  the  stress  level. 
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4 . 6   Poisson's  Ratio 

Although  Poisson's  ratio  v  and  modulus  of  elasticity 
E  are  both  important  for  calculations  of  stresses  and  deform- 
ations in  soil,  the  effect  of  variations  in  Poisson's  ratio 
has  received  very  little  attention  and  this  parameter  has 
frequently  been  assumed  to  be  constant.   For  a  uniaxial 
stress  increment,  Poisson's  ratio  is  defined  as  the  ratio  of 
lateral  to  axial  strain. 

£3 
v  =   -   -i  (4.25) 

el 

in  which  v  =  Poisson's  ratio,  £_.  and  e,  are  the  lateral  and 

axial  strains  respectively.   Chen  (1948)  found  that  the 

value  of  Poisson's  ratio  for  cohesionless  soils  varies  with 

strain.   He  defined  the  tangent  Poisson's  ratio  as 

Ae_ 
vt=   -   ^1  (4.26) 

where  v  denotes  tangent  or  incremental  Poisson's  ratio, 
Ac,  and  Ae,  are  incremental  values  of  e-.  and  e...  Tangent 
Poisson's  ratio  as  given  by  equation  (4.2  6)  has  been  used 
by  several  authors.  Duncan  and  Chang  (1970)  used  the  same 
equation  in  a  different  form  for  incremental  analysis,  in 
which  incremental  Poisson's  ratio  has  been  defined  as 


v,_  = 


Ae,  -Ae 
1    v 


t  "     2Ae,  (4.27) 


in  which  Ae   is  the  increment  in  volumetric  strain.   They 
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found  values  of  Poisson's  ratio  ranging  from  0.11  to  0.65. 
Values  of  tangent  Poisson's  ratio  may  range  from  as  low  as 
0.1  to  as  high  as  1.6. 

Poisson's  ratio  can  be  expressed  in  terms  of  shear 
modulus   G  and  bulk  modulus  K,  which  can  be  written  as: 


v  = 


3K  -  2G 
6K  +  2G 


(4.28) 


In  some  cases  (vibration  analysis)  it  is  more  convenient  to 
obtain  K  and  G,  whence  v  can  be  calculated  from  equation 
(4.28) 

Lade  (1972)  presented  a  formulation  of  tangent 
Poisson's  ratio  assuming  the  hyperbolic  stress-strain  func- 
tion proposed  by  Duncan  and  Chang   (1970) . 


G  -  Flog( — -) 


Vt   = 


(a. 


1  - 


a3) 


o  ~   n 

K-p  ( -) 

a      P= 


1  - 


Rf  •  (c^-c^)  (1-SincJ)) 
2  OCoscj)  +  2  a3«Sin<j) 


(4.29) 


in  which  G  is  the  value  of  tc.ngent  Poisson's  ratio  v.  at 
o-  =  1  atmosphere,  F  is  the  decrease  in  value  of  v.  per  log- 
cycle  increase  in  o,  and  d  is  a  constant  of  hyperbolic  vari- 
ation,  other  parameters  have  been  defined  in  equation 
(4.4).   Lade  (1972)  showed  that  a  plot  of  -  e3  and  c1   for  a 
given  confining  pressure  on  a  log-log  diagram  approximates  a 
straight  line  (Figure  4.6a).   This  straight  line  can  be 
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(Data  From    Duncan   and    Chang,  1970) 
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expressed  as  an  expotential  function 

-  e3  =  P  •  e^  (4.30) 

where  P  is  the  value  of  -  £-.  at  e,  =  1%  and  m  is  the  slope 
of  the  straight  line.   For  some  soils  similar  plots  for  var- 
ious confining  pressures  show  more  or  less  the  same  slope 
of  the  line;  in  orther  words  value  of  m  is  independent  of 

confining  pressure.   The  position  of  the  line  is  shifted 

with  confining  pressure;  that  is,  values  of  P  will  be  dif- 

°3 
ferent.   If  values  of  P  are  plotted  against  ( — )  in  a  log-log 

pa 

diagram,  it  results  in  a  straight  line  (Figure  4.6b).   This 

line  can  be  expressed  as: 

a  ■»  q 

P  =  L  (—  )  (4.31) 

Pa 

°3 
in  which  L  is  the  intercept  at  —  =  1  and  q  is  the  slope 

"a 

of  the  line.   Substituting  for  P  from  equation  (4.31)  into 

equation  (4.30)  results 

a*     q   m 

-  e,  »  L  (—  )    e™  (4.32) 

Tangent  Poisson's  ratio  can  be  determined  from  the  above 
equation  (4.32)  by  taking  the  first  derivative  with  respect 
to  £..  . 


de_  a-,   q        , 

v.    =   -  -r-3-  =  L  •  m  •(  -^  )    •£?  (4.33) 

tan       de1  Pa         1 
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value  of  e1  in  this  equation  is  the  total  major  principal 
strain  in  percent.   Using  parameters  L,  m,  q  and  values  of 
e,    the  tangent  Poisson's  ratio  for  primary  loading  can  be 
determined  for  any  state  of  stress  in  triaxial  compression. 
Lade  (1972)  has  shown  good  agreement  of  the  experimental  and 
predicted  values  of  tangent  Poisson's  ratio  up  to  about  2% 
axial  strain. 

From  the  results  of  triaxial  tests  the  coefficients 
in  equations  (4.4)  and  4.33)  can  be  evaluated  and  values  of 
E   and  v   can  be  obtained  for  plane  stress  conditions.   In 
principle,  the  corresponding  values  for  plane  strain  condi- 
tions should  be  given  by 

Et 
E    =   ^_  (4.34) 

pl     (1  -vj) 

v   =  — ^ 
Pl   '  l-V-t  (4.35) 

but  the  carefully  conducted  triaxial  and  plane  strain  tests 
by  Lade  (1972)  showed  that  this  is  not  the  case  (Figures  4.6c 
and  4.6d) .   For  this  reason,  it  was  decided  to  fit  spline 
functions  to  the  actual  test  data  and  obtain  E   and  v, 
directly  from  the  data  for  any  given  stress  state,  which  is 
defined  in  terms  of  octahedral  stresses.   Details  of  the 
procedure  are  given  in  Section  4.9.   It  is  recognized  that  an 
incrementally  elastic  representation  of  soil  properties  is  in- 
accurate if  the  soil  is  approaching  failure,  or  is  in  a  post- 
peak  stress-strain  range,  because  for  these  conditions  plastic 
strains  dominate  behavior.   However,  suitable  elastic-plastic 
formulations  are  still  in  too  formative  a  stage  (Lade  and 
Duncan,  1976)  for  practical  use. 
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4. 7   Anisotropy 

Strength  anisotropy  in  vertical  and  horizontal  direc- 
tions of  naturally  deposited  soil  is  well  accepted  in  the  lit- 
erature.  The  reason  for  strength  anisotropy  is  that  during 
deposition  of  a  soil  layer  by  a  natural  process,  soil  par- 
ticles attain  some  preferred  orientation.   It  is  also  known 
that  soil  particles  do  not  have  any  regular  shape.   As  a 
result,  a  soil  layer  offers  more  resistance  to  deformation 
due  to  an  imposed  load  in  one  direction  than  another  direc- 
tion.  Madhav  and  Roy  (1970)  have  determined  the  anisotropy 
ratio  of  vertical  to  horizontal  shear  strength  about  1.5. 
Holubec  (1968)  has  shown  that  for  cohesionless  soil  Young's 
modulus  in  vertical  direction,  E   is  greater  than  that  in 
horizontal  direction,  E„  ,  i.e.,  E   >  E„  as  shearing  in- 

n  V      n 

creases.   For  initial  equal  all  around  pressure  Poisson's 
ratio  in  vertical  direction,  v   is  greater  than  that  in 
horizontal  direction  v„  .   ie.,  vtt  >  v„  ;  and  as  shearing 

n  V     n 

increases,  a  reversal  occurs,  vtt  <  v„.   If  E  >E  ,  no  matter 

V      ri  V   rl 

what  is  the  relation  between  v  .  and  v„  ,  horizontal  stress 

V        n 

will  be  less  than  if  E   =  E„. 

V      n 

Poulos  (1972)  has  shown  that  horizontal  deflection  of 
a  foundation  can  not  be  predicted  assuming  soil  as  homo- 
geneous, elastic  isotropic  material.   He  cited  the  following 
parameters  as  having  important  effects  on  the  mechanism  of 
stress  and  strain  distribution  in  a  soil  medium:   (1)  Hori- 
zontal soil  modulus  value,  E„  and  the  ratio  of  E  /E  , 
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(2)   Horizontal  Poisson's  ratio,  \>  has  greater  influence, 

2 
as  vertical  deflection  is  a  function  of  (1-v  )  whereas 

horizontal  deflection  is  a  function  of  (1+v)  (l-2v) . 

Wright  and  Duncan  (1972)  studied  sliding  of  the  Waco 
Dam.   Investigation  showed  that  undrained  shear  strength  of 
the  material  was  highly  anisotropic.   Ratio  of  horizontal 
to  vertical  strength  was  about  0.4.   Analyses  based  on  the 
assumption  that  the  soil  is  isotropic  indicates  a  stable 
condition  of  the  dam,  whereas  use  of  anisotropic  strength  as 
determined  by  experiments,  produced  results  in  agreement  with 
the  observed  failure.   This  study  shows  the  significance  of 
considering  anisotropy  in  analysis. 

The  most  commonly  used  pipe  materials  are  concrete, 
steel  and  aluminum,  which  can  reasonably  be  considered  as 
isotropic  materials  in  the  stress  range  preceding  failure. 
Hence  no  consideration  of  anisotropy  is  given  to  these 
materials.   For  pipes  constructed  of  composite  plastic  lam- 
inates, the  effects  of  anisotropic  properties  are  important, 
and  can  be  handled  by  following  the  procedure  described  in 
the  next  section. 

Anisotropic  Pipe  Properties  (for  thick  pipe  only) 
Pipes  constructed  of  composite  plastics  may  be 
highly  anisotropic;  also,  the  directions  of  principal 
anisotropy  change  continually  along  the  perimeter  of  pipe. 
This  requires  transformation  of  elastic  properties  according 
to  its  geometric  orientation. 


88 


Stavsky  and  Hoff  (1969)  presented  a  procedure  for 
transformation  of  elastic  properties  for  orthotropic  compos- 
ite materials.   The  method  is  based  on  the  theory  of  elas- 
ticity and  was  presented  in  detail  by  Lakhnitskii  (translated 
from  Russian  by  P.  Fern,  1963).   In  Figure  4.7  a  block  of  an 
orthotropic  material  is  shown  along  with  the  global  coordi- 
nate direction  (x,y) .   At  any  point  the  normal  direction  N, 
and  the  tangential  direction  S,  is  shown  at  an  angle  a  with 
the  global  coordinate  axes.   In  the  following  developments 
subscripts  x  and  y  will  denote  quantities  in  X  and  Y  direc- 
tions.  Similarly  subscripts  n  and  s  will  be  used  to  refer 
to  quantities  in  normal,  N  and  tangential,  S  directions. 

According  to  the  theory  of  elasticity,  in  the  most 
general  case  the  expression  for  elastic  potential,  V,  can  be 
considered  in  terms    of   component  stresses  and  strains, 
and  can  be  written  as: 

V  =  T(a  e   +a  e   +  ...  +  t   y   )  (4.36) 

2   x  x   y  y  xy '  xy 

in  which  a  and  e  refer  to  normal  stress  and  strain,  x      and 
Y   represent  shear  stress  and  shear  strain.   Then  the  elastic 
potential  is  expressed  in  terms  of  rotated  coordinate  system. 
As  the  elastic  potential  is  independent  of  the  coordinate 
system,  the  two  expressions  can  be  equated.   Then,  equating 
the  coefficients  of  similar  terms,  the  unknown  constants  of 
transformation  can  be  determined.   After  simplification,  the 
following  set  of  equations  can  be  written,  which  relate 
properties  in  the  transformed  axes  to  those  in  the  global  axes. 
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FIGURE   4.7     COORDINATE    SYSTEM    FOR    TRANSFOR 
MATION    OF  ANISOTROPIC    MATERIAL 
PROPERTIES. 
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Sin   2a 


(4.  37a) 


—     .  4     En  4     1 

E   =  Sin  a   +   =7^  •  Cos  a  +  -r 

x            E  4 
s 


n     ~ 

rT-       ~     2v 
G         nS 
ns      " 


Sin   2a 


(4.  37b) 


Vyx  =  E* 
J      n 


K,s  -  i(1  +  2  -ns  +  f1  -  !r->  •  sin2  2a 

s    ns 


(437c) 
in  which  G  and  v   denotes  shear  modulus  and  Poisson's 
ratio- 
According  to  Maxwell's  reciprocal  theorem  the  two 

strains  e   and  e  must  be  equal  if  the  two  stresses  a     and 
n      s  n 

o   are  of  equal  magnitude  and  sense.   Consequently, 


V  E 

ns  .  _n 

v  E 

sn  s 


(4.37d) 


Shear  strain  v    due  to  stress  a   is 
'yx  y 


'  mlay 
"yx     E 


(4.37e) 


n 


in  which  m.  is  given  by 


m,  =  Sin  2a  • 


E      E 
n     n 


^ns  '  Es    2'Gns 


2 
Cos  a 


(1  +  2v 


ns 


E     E 

E     G    ' 
s    ns 


(4.37f ) 
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Shear  stress  t    will  produce  a  component  of  strain 


e  ,  which  is  given  as 


m„x 

= 2  y* 


(4.37g) 


where  m_  is  given  by 


m~  =  Sin  2a 


Vns  +  E 


n 


n 


2G 


-  Sin  a  • (1  +  2v 


ns 


ns 


n 


n 


ns 


(4.37h) 


In  the  above  equations  G    is  given  as 
^         ns     ^ 


ns 


yx 


ns 


n 


(  Y   +    2vns  +  IT  ) 
s 


(  1  +  2v    +  . 
ns    E 


n 


E 


n 


ns 


)  -  Cos   2a 


(4.  37i) 


In  equation  (4.37a)  through  (4.37i)   E  ,  E   and  v 

will  be  defined  for  an  orthotropic  material.   Values  of 

E  ,  E   G   will  be  calculated  for  a  given  element  depending 
x'   y   yx  3 

on  its  position  on  the  pipe. 

The  next  step  is  to  write  the  equations  for  the  gen- 
eralized Hook's  law  for  plane  strain  conditions  in  the  global 
coordinate  system.   In  matrix  form  these  equations  can  be 
written  as 
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'                        ' 
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Yxy 

1        yx 
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Y 
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m„  _, 


n 


m. 


n 


yx 


i  yx    J 


(4.38  ) 


or 


{a}  =  [D]  {e} 


(4.39) 


Inversion  of  the  above  equation  (4.39  )  will  result 
in  the  appropriate  [D]  matrix,  which  will  be  used  as  the 
elasticity  matrix  in  equation  (3.14a)  to  generate  the 
stiffness  matrix  of  the  orthotropic  pipe  element. 

4 . 8   Determination  of  Elasticity  Matrix  [D]  for  Soil 

The  stiffness  matrix  for  a  finite  element  consists 
of  two  key  ingredients,  one  of  them  is  geometry  and  dis- 
placement function,  the  other  is  material  property  and  type 
of  analysis.   In  equation  (3.14)  [b]  matrix  takes  care  of 
geometry  and  shape  function;  [D]  matrix  reflects  the  material 
properties  and  type  of  analysis.   In  this  section  development 
of  [D]  matrix  considering  anisotropy  and  plane  strain  condi- 
tions is  presented. 

The  general  stress  strain  relations  for  an  elastic 
anisotropic  medium  can  be  written  by  assuming  material 
properties  are  known  in  two  mutually  perpendicular  directions 
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(namely  vertical,  y  and  horizontal  x  and  z  directions). 

(Lekhnitskii,  1963) 

a  a  a 

x  y  z 

Ex  "  E    "  Vyx*E    "  Vxy*E  (4.40a) 

x      J         y  *      x 

°x    av  CTz 

e   =  -  v   —  +  =*~  v  •  ~               (4.40b) 

y       yx  E     E      yx  E 

J       J    x     y     J  x 

a        a    a 

e   =  -  v   •=—  -  v   »=^  +  =—  (4.40c) 

z       xy  E      yx  E     E 

J  x       y    x 

Txz 
Y    =   2(1  +  v   )  — ~  (4.40d) 

'xz  xy    E 

J     x 

yxy  =  2(1  +  V"i£  {4-40e) 

V  =  2(1  +  Vir  (4-40f) 


For  plane,  strain  condition  e   =0  and  also  v    an<3 
r  z  '  xz 

v    are  zero.   a   can  be  determined  in  terms  of  a   and  a 
'yz  z  x      y 

from  equation  (4.40c)  for  e   =0.  a      in  equations  (4.40a) 

and  (4.40b)  can  be  eliminated.   Ratio  of  moduli  can  be  de- 

E  G 

fined  as  modular  ratios  =^  =  r  and  m  =  =£%-     where 

is  hi 

y  y 

E   E 

_* y_ 


Jxy    E   +  E     2v   E 
1  x    y  +    yx  x 

Using  these  definitions,  the  stress-strain  relations 
can  be  written  in  matrix  form  similar  to  equation  (4.38). 
The  [D]  matrix  can  be  determined  by  inverting  the  stress- 
strain  matrix  and  can  be  written  as: 
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[D]  = 


where 


(1+v   )  •  (1-v    -  2-rv   ) 
*    xy      xy       yx 


b 

0 


0 
0 
d 


(4.41) 


a  =  r  •  (1  -  rv    ) 

yx 


b  =  r  •  v    (1  +  v    ) 
yx       xy 


c  =  (1  -  v   ) 
xy 


d  =  m  •  (1  +  v   )«(l-v   -  2r  •  v   ) 
xy        xy         yx 

The  expression  for  [D]  as  given  by  equation  (4.41) 

considers  both  anisotropic  soil  properties  and  plane  strain 

condition.   In  the  above  equation  E  will  be  determined 

using  equation  (4.18);  and  corresponding  Poisson  '  s  ratio  v   ; 

xy 

v    from  (4.37d)  and  values  of  r  will  be  defined  as  input 
yx 

data.   This  completes  the  evaluation  of  [D]  matrix,  which 
will  be  used  in  equation  (3.14)  to  evaluate  stiffness 
matrix  for  soil  elements. 


4 . 9   Use  of  Spline  Function  for  Representation  of  Material 
Properties 

Use  of  spline  function   for   representing   non- 
linear, stress  dependent  properties  for  soil  have  been  pre- 
sented in  sections  4.5  and  4.6.   In  this  section  a  procedure 
for  reduction  of  test  data  is  presented,  which  requires  no 
assumptions — such  as  parabolic,  hyperbolic,  or  exponential-- 
to  obtain  a  "fit"  to  the  test  data. 


95 


Lade  (1972)  conducted  plane  strain  tests  on  loose 

sand  (Monterey  No.  0)  for  several  values  of  confining 

pressure  and  presented  results  of  deviator  stress  (a-,  -aJ  , 

major  principal  strain,  e,  and  volume  strain,  e  .  Figure 

4.8  shows  the  stress-strain  and  volume  change  curves  that 

were  obtained.   In  this  study  octahedral  stresses  are  used, 

so  it  is  necessary  to  change  the  test  data  from  principal 

stress  to  an  octahedral  stress  system.   Every  point  on  a 

curve  has  unique  value  of  a,  and  a^.   Value  of  o.  can  be 

determined  for  plane  strain  conditions  from  equation  (4.40) 

as  values  of  e,  can  be  determined  from  e,  and  e      and 

e_  =  0.    For,  each  point,  octahedral  normal  stress  a   . , 

and  octahedral  shear  stress  t   .  are  obtained  from  equations 

oct 

(D.l)  and  (D.2),  and  equation  (D.8)  can  be  used  to  evaluate 

t   L/,  at  failure  for  a  given  a   ..   The  stress  ratio  at 
oct/f  oct 

OCt/ 


failure  t   .  /       is  determined. 
-A 


oct/f 


Procedure  for  Data  Reduction 

(1)  Select  successive  values  of  (o^   -a3) /  ev  and  £^ 
from  the  test  data  at  a  given  o^   starting  from  zero. 

(2)  Cubic  splines  are  fitted  to  the  data  following 
the  procedure  described  in  section  4.5. 

(3)  Values  of  tangent  modulus  E  ,  and  tangent 
Poisson's  ratio,  v  ,  are  evaluated  using  equation  (E.16) 
and  (E.10).   (See  Appendix  E) 
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cr  =  1.2  kg  /cm 


0.4 


FIGURE  4.8     PLANE    STRAIN    TEST  ON    LOOSE 
MONTEREY    NO.  O  SAND. 
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(4)  Small  increments  of  e,  are  chosen  and  for  each 

e, ,  values  of  a, ,  a„,  o,,  a   . ,    x      .,    t   .  ,_    ,  ..   ^  . , 
1  1    2'      3    oct   oct   oct/f  and  the  failure 

stress  ratio  are  evaluated  using  procedures  described  above 

and  in  section  4.5. 

(5)  Step  3  is  repeated  for  each  increment  of  e,  up 
to  failure. 

(6)  Step  1  through  5  are  repeated  for  test  data  at 
different  values  of  a... 

(7)  All  values  of  E.  are  plotted  against  their 

corresponding  values  of  a  and  values  of  stress  ratio 

r     3  oct 

t   ,  .       are  noted  for  all  points. 

OCt/T    ....  ^ 

'    oct/f 

(8)  Contours  are  drawn  for  selected  values  of  the 

stress  ratio.   For  convenience,  values  of  stress  ratio  should 

range  from  zero  to  unity  with  increments  of  0.1.   This  will 

generate  11  curves  of  E.  vs.  a  for  stress  ratio  ranging 

'  t      oct  3 

from  0.0  to  1.0  (Figure  4.9). 

(9)  Steps  7  and  8  are  repeated  for  tangent  Poisson's 
ratio  v   (figure  4.10). 

(10)   From  the  peak  points  in  Figure  4.8  values  of 

t   ,_  ,,  and  o   ,_  ...  are  calculated  and  shown  plotted  in  Figure 
oct/f      oct/f  r 

4.11. 


Lade  (1972)  also  conducted  plane  strain  tests  on 
dense  sand.   Figure  4.12  shows  typical  stress-strain  and  vol- 
ume change  results.  Following  the  procedure  outlined  above 

values  of  E,_  and  v^  are  plotted  against  a     for  different 
t      t     *  oct 
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FIGURE   4.9     TANGENT    MODULUS    VS.   OCTAHEADRAL 
NORMAL    STRESS    AND   FAILURE    RATIO 
FOR    LOOSE    SAND. 
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FIGURE  4.10  TANGENT  POISSON'S  RATIO  VS.  OCTA- 
HEADRAL  NORMAL  STRESS  AND  FAI- 
LURE   RATIO   FOR    LOOSE   SAND. 


100 


CM 


I I I I I I I 


I I L 


UJD/6)|    *43°J. 


O 


< 

Q 

< 
tu 

iS 

o 
o 

>§ 

ecu 

wo 
o 

a:-1 

<  « 

ujuj 

x£ 

^2 

ho: 

yp 

Ocn. 


id 

Z> 


101 


0.4L 


FIGURE   4.12   PLANE    STRAIN    TEST   ON    DENSE 
MONTEREY    NO.  0    SAND.  (After 
Lade ,  1972 ) 
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values  of  stress  ratio  at  failure  and  shown  in  Figure  4.13 

and  4.14,  respectively.   Figure  4.15  shows  the  plot  of 

t   ,  ,j-   and  a  .  ,c    for  dense  sand, 
oct/f      oct/f 

If,  as  in  the  plane  strain  test,  the  principal  planes 
remain  fixed  then  the  procedure  described  previously  is  a 
precise  general  method  for  representing  the  constitutive 
relations  for  granular  soils  and  storing  the  results  in  the 
computer.   Only  states  of  stress  below  failure  are  consid- 
ered, and  loading  must  increase  continuously.   However,  if 
the  principal  planes  rotate  during  loading,  the  directions 
of  the  principal  stress  and  principal  strain  increments  may 
not  coincide.   The  potential  effects  on  the  stress-strain 
relation  is  not  accounted  for  in  this  study. 
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(Tocf  kg  /cm 


FIGURE  4.13  TANGENT  MODULUS  VS.  OCTAHEADRAL 
NORMAL  STRESS  AND  FAILURE  RATIO 
FOR    DENSE    SAND. 
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CHAPTER  V 
COMPUTER  PROGRAM 

5 . 1   Introduction 

In  Chapter  II  the  need  for  a  realistic  analytical 
tool  for  the  solution  of  culvert  problems  was  described  in 
some  detail.   The  vehicle  for  solution,  that  is,  the  finite 
element  method  was  discussed  in  Chapter  III  and  methods  for 
considering  representative  material  properties  were  treated 
in  Chapter  IV.   It  is  apparent  that  the  number  of  variables 
encountered  in  culvert  problems  can  be  handled  only  with  the 
help  of  a  large  digital  computer,  and  that  an  efficient  tech- 
nique for  numerical  analysis  as  required.  In  this  chapter, 
computer  program,  data  structure,  and  solution  technique 
developed  for  this  study  will  be  presented  with  the  aid  of 
simplified  flow  charts. 

The  Computing  Center  at  Purdue  University  has  a  CDC 
6500/6400  system  computer.   Therefore,  the  program  has  been 
written  using  features  of  CDC  system.   The  language  is 
FORTRAN  IV,  and  the  program  is  called  FINLIN,  (Finite 
element,  Isoparametric,  Non-Linear  with  Interaction  and  No- 
tension)  . 
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The  following  attractive  features  were  utilized  in 
developing  the  computer  program: 

(1)  Isoparametric  triangular  finite  elements  with 
one  curved  boundary,  and  independent  of  the  local  coordinate 
system,  was  used  to  represent  soil  materials.   The  displace- 
ment function  is  quadratic  resulting  in  linear  strain  dis- 
tribution within  the  element. 

(2)  Nonlinear,  stress  dependent  anisotropic  soil 
properties  can  be  used.   Tangent  Young's  modulus  and  tangent 
Poisson's  ratio  for  any  state  of  stress  is  determined  by 
fitting  spline  functions  to  actual  test  data  which  are 
represented  in  terms  of  octahedral  stresses. 

(3)  A  soil-pipe  interaction  element  was  introduced 
consisting  of  a  curvilinear  rectangle  of  zero  thickness, 
i.e.,  two  adjacent  corners  of  the  rectangle  have  the  same 
coordinates. 

(4)  Bilinear  stress-strain  behavior,  which  includes 
slip  at  a  critical  ratio  of  shear  stress  to  normal  stress, 
was  used  to  represent  soil-pipe  interface  behavior. 

(5)  A  "no-tension"  analysis  was  used  both  for  soil 
and  interaction  elements. 

(6)  Realistic  representation  of  pipe,  including 
rotational  stiffness  at  nodes  has  been  introduced  using  a 
segment  of  a  curved  bar  with  six  degrees  of  freedom. 

(7)  Self-weight  of  soil  is  distributed  throughout 
the  soil  mass,  and  different  types  of  materials  in  different 
zones  can  be  handled. 
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(8)  By  selecting  appropriate  options  analysis  using 
linear  or  nonlinear  soil  properties  can  be  performed. 

(9)  The  program  is  capable  of  incremental  analysis 
in  layers,  which  realistically  models  the  construction  se- 
quence. 

Limitations  of  the  computer  program  are  as  follows: 

(1)  Only  plane-strain  conditions  have  been  con- 
sidered, so  that  three-dimensional  effects  are  ignored. 
However,  the  program  can  be  used  for  two-dimensional  prob- 
lems in  plane  stress  conditions. 

(2)  Time-dependent  soil  properties  and  prestressing 
due  to  compaction  has  not  been  accounted  for. 

(3)  Yielding  or  cracking  of  the  pipe  has  not  been 
considered. 

(4)  Only  static  loads  can  be  analyzed. 

(5)  There  is  no  provision  for  automatic  mesh  gener- 
ation. 

5. 2   Program  Outline 

Any  finite  element  scheme  requires  the  following  se- 
quence of  operations. 

(1)  Reading  information  about  nodes,  elements  and 
boundary  conditions. 

(2)  Input  of  material  properties. 

(3)  Generation  of  element  stiffness  matrix  and 
assembly  to  form  the  structure  stiffness  matrix  with  modifi- 
cations for  given  boundary  conditions. 
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(4)  Solution  of  stiffness  equations  for  given  inter- 
nal and  external  loads. 

(5)  Determination  of  element  strains  and  stresses. 

(6)  Recycling  according  to  given  options  for 
analysis. 

The  computer  program  FINLIN  follows  a  similar  se- 
quence of  operations.   Descriptions  of  each  of  these  steps 
are  presented  in  later  section.   First,  some  special  fea- 
tures of  the  computer  system  will  be  discussed,  which  are 
helpful  in  writing  programs  which  require  large  memory  loca- 
tions and  involve  very  large  numbers  of  numerical  operations. 

Overlay 

Overlay  is  a  high  speed  loader  system  which  trans- 
fers relocatable  program  sources  to  central  memory  for  exe- 
cution.  Each. overlay  consists  of  a  main  program  and  associ- 
ated subroutines.   Each  overlay  is  identified  with  a  pair  of 
intiger  numbers,  the  first  number  is  for  primary  and  the 
second  for  secondary  levels.   Transfer  of  all  information 
from  one  overlay  to  another  can  be  achieved  by  using  blank 
common,  labelled  common  or  by  disc  tape  units.   All  overlays 
excepting  the  main  overlay  (0,0)  are  stored  outside  the  cen- 
tral memory.   Any  primary  or  secondary  overlay  can  be  loaded 
to  the  memory  when  called  by  an  appropriate  control  card. 
The  main  advantage  of  using  an  overlaid  system  is  that  por- 
tions of  a  program  which  are  not  needed  at  one  time  can  re- 
main outside  the  central  memory  and  does  not  occupy  these 
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valuable  spaces.   As  a  result,  more  storage  space  is  avail- 
able for  usage.   For  example,  when  the  program  is  evaluating 
the  stiffness  matrix  and  solving  for  displacements,  routines 
for  stress  or  strain  are  not  needed.   Hence,  routines  for 
stress  and  strain  can  be  stored  elsewhere. 

(Disc)  Tape  Unit 

These  units  are  used  for  temporary  READ/WRITE 
operations  for  unformated  and/or  formated  operations.   Large 
amounts  of  information  (e.g.,  nodal  and  element  data),  which 
are  used  several  times  in  the  program  can  be  written  on  tape 
units  and  the  memory  locations  thus  released  can  be  used  for 
some  other  purpose.   Necessary  sets  of  data  can  be  read  from 
the  tape  unit  when  required;  a  desired  set  of  data  can  be 
found  by  following  the  sequence  in  which  the  data  has  been 
written  on  the- unit.   So  a  bookkeeping  of  the  sequence  has 
to  be  maintained  throughout.   Also  READ  and  WRITE  operations 
can  not  be  mixed.   Though  rewinding  a  tape  unit  is  permitted 
at  any  stage,  the  READ  or  WRITE  operations  have  to  be  se- 
quential.  This  limitation  causes  some  difficulty  for  general 
use. 

Random  Mass  Storage 

This  technique  is  similar  to  the  disc  tape  unit  but 
the  need  for  bookkeeping  has  been  eliminated.   This  is  done 
by  assigning  an  index  number  for  each  record  and  defining 
the  length  of  each  record.   The  READ/WRITE  operations  can  be 
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arbitrary  rather  than  in  sequential  form.   This  facility  can 
be  used  very  efficiently  for  nonlinear  analysis  where  the 
state  of  stress  and  strain  are  required  very  often.   The 
element  numbers  can  be  used  as  index  number. 

The  above  mentioned  features  have  been  fully  utilized 
in  writing  the  computer  program.   The  scheme  of  operations 
has  been  subdivided  into  several  OVERLAY'S  whose  functions 
are  described  here. 

Figure  5.1  shows  the  organization  of  a  complete  pro- 
gram which  consists  of  several  OVERLAYS.   Figure  5.2  shows 
the  flow  chart  for  OVERLAY  (0.0)  which  happens  to  be  the 
main  OVERLAY.   Calling  and  execution  of  all  other  OVERLAYS 
are  controlled  by  this  section  depending  upon  the  type  of 
analysis  desired. 

The  flow  diagram  for  OVERLAY  (1,0)  is  shown  in 
Figure  5.3.   Description  of  problem  geometry,  materials, 
number  of  nodes  and  elements  are  read  here.   Semi-band 
width  of  the  stiffness  matrix,  number  of  degrees  of  freedom 
are  calculated.   Nonlinear  material  properties  are  read  for 
each  type  of  material.   Data  reduction  and  spline  fitting  is 
performed  and  pertinent  parameters,  such  as  tangent  modulus 
and  tangent  Poisson's  ratio,  are  evaluated  and  stored  for 
future  use.   Geometry,  properties  and  location  of  pipe  is 
also  read  here.   Evaluation  of  area  for  triangular  elements, 
and  other  values  are  computed  here. 

OVERLAY  (2.0)  plots  the  finite  element  mesh  from  the 
data  read  earlier,  Figure  5.4. 
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OVERLAY  (6,0) 
MODIFY  STIFFNESS  MATRIX  FOR  BOUNDARY 
CONDITION  AND  SOLVE  FOR  INCREMENTAL 

DISPLACEMENTS 


Figure  5.1  General  Flow  Chart  of  program  FINLIN 
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CALL  OVERLAY  (3,0) 


CALL  OVERLAY  (6.0) 
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CALL  OVERLAY  (4.0) 
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CALL  OVERLAY  (5.0) 


<END  OF  NONLINEAR  ANALYSIS?> 
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Figure  5.2.  Flow  Chart  for  OVERLAY  (0,0) 
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READ  NO.  OF  NODES,  ELEMENTS,  MATERIALS, 
LOADS  AND  CONTROL  PARAMETERS 


READ  NODAL  AND  ELEMENT  DATA 
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READ  PIPE  PROPERTIES,  GEOMETRY,  LOCATION 


CALCULATE  SEMI-BAND  WIDTH,  CALCULATE 
ELEMENT  INFORMATION 


READ  AND  STORE  EXTERNAL  LOADS 


PRINT  INPUT  INFORMATION 


READ  STRESS-STRAIN  AND  CONFINING 
PRESSURE  DATA,  TYPE  OF  TEST 


FIT  SPLINE  FUNCTION,  FIND  SPLINE 
PARAMETERS 


DETERMINE  TANGENT  MODULUS  AND  POISSON'S 
RATIO  AT  PRESET   STRESS  RATIO 
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CONSIDERED?, 
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(continue) 


Figure  5.3   Flow  Chart  for  OVERLAY (  1,0  )  ,  (  Continued  ) 
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Figure  5.3   Flow  Chart  for  OVERLAY  (1,0) 
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USE  CALCOMP  ROUTINES  FOR 
PLOTTING  FINITE  ELEMENT  MESH 


Figure  5.4  Flow  Chart  for  OVERLAY  (2,0) 
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The  most  important  portion  of  the  program  is  OVERLAY 
(3.0),  which  is  shown  in  Figure  5.5.   In  this  section  evalu- 
ation of  the  stiffness  matrix,  structural  assembly,  and  esti- 
mation of  nodal  loads  is  carried  out.   Because  of  limited 
capacity  of  the  computer's  core  memory,  all  elements  in  the 
structure  can  not  be  assembled  at  one  time.   This  is  done  in 
blocks  the  size  of  which  should  be  at  least  a  square  array 
of  the  size  of  semi-band-width  of  the  assembly  [semi-band- 
width is  evaluated  and  stored  in  OVERLAY  (1,0)].   Then  a 
searching  scheme  locates  which  of  the  nodes  and  corresponding 
elements  should  fit  in  the  first  block.   Evaluation  of 
stiffness  matrix  is  performed  element-by-element.   The 
three  types  of  elements,  curved  bar,  interaction,  and  tri- 
angular elements,  are  each  treated  separately  because  of 
their  different  characteristics,  and  are  described  here. 

Type  1,  Curved  Bar  Element 

The  pipe  radius,  EI  value,  angle  B,  (Figure  B. 1) , 
position  of  center  of  pipe  and  angles  <}>,  and  c(>2  (Figure  B.2) 
need  to  be  calculated  from  the  information  supplied  by 
OVERLAY  (1,0).   As  properties  of  the  pipe  remain  constant, 
no  particular  attention  is  required  for  non-linear  behavior. 
Stiffness  coefficients  for  element  1-2  (Figure  B.l)  are 
evaluated  and  transformed  to  the  global  axes  using 
equation  (B. 16) . 
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SET  BLOCK  SIZE  FOR  STRUCTURAL 
STIFFNESS  MATRIX  AND  LOAD  VECTOR 


FIND  NUMBER  OF  ELEMENTS  IN  BLOCK 


READ  STRESSES  AND  STRAINS  FROM 
PREVIOUS  CYCLE,  CHOOSE  NUMBER 
AND  INTERVAL  OF  SPLINE 


FIND  STRESS  DEPENDENT,  ANISOTROPIC 
PROPERTY,  APPROPRIATE  POISSON'S 

RATIO 


FIND  APPROPRIATE  STIFFNESS  MATRIX 


JO INT -ELEMENT  STIFFNESS  TO  FORM 
STRUCTURAL  STIFFNESS  MATRIX 


STORE  STIFFNESS  MATRIX  AND  LOAD 
VECTOR  ON  MASS  STORAGE  UNIT 


IF  ALL 
ELEMENTS 
CONSIDERED? 


Yes 


Figure  5.5.   Flow  Chart  for  OVERLAY  (3,0) 
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Type  2,  Interaction  Element 

There  are  two  alternatives  depending  whether  inter- 
action behavior  is  desired  or  not.   If  interaction  is  not 
required,  properties  of  interaction  are  maintained  constant 
throughout  and  evaluation  of  stiffness  matrix  proceeds  as 
described  in  Appendix  C. 

If  interaction  behavior  is  required  then  for  the 
first  step  the  scheme  follows  the  procedure  for  a 
no-interaction  situation.   In  subsequent  steps,  the  normal 
stress  on  the  element  due  to  all  previous  increments  are 
read  from  mass  storage  data.   If  the  normal  stress  is  com- 
pressive, and  the  ratio  of  shear  stress  to  normal  stress 
is  less  than  the  limiting  value,  the  modulus  in  tangential 
direction,  E   is  maintained  at  the  same  high  value  as  that 

for  normal  direction,  E 

n 

E   =  E   +  °°  (5.1) 

n    t 

If  the  normal  stress  on  the  interaction  element  is 

zero   in   nature,  both  E^_  and  E   are  set  to  small  values 

t      n 

(because  zero  values  for  E   and  E^.  will  create  ill- 

n      t 

conditions) . 

If  the  normal  stress  on  the  interaction  element  is 
compressive  but  the  ratio  of  shear  stress  to  normal  stress 
is  larger  than  the  limiting  value,  the  E   value  is  kept  un- 
changed while  E.  is  set  to  a  small  value,  which  simulates 


slip. 
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In  case  of  failure  of  the  interaction  element  either 
due  to  excessive  shear  stress  or  to  tension,  a  message  will 
be  printed  to  that  effect.   In  case  of  tension,  excess 
tension  is  eliminated  by  "no-tension"  analysis,  as  described 
later. 

Triangular  Element 

This  type  of  element  is  used  for  the  soil  mass.   It 
can  have  linear  elastic  or  nonlinear  properties  depending 
upon  the  option.   Anisotropy  parallel  to  the  general  coord- 
inate system  is  also  permitted  and  can  be  handled  as  de- 
scribed in  section  4.8.   For  linear  elastic  properties,  the 
modulus  value  in  vertical  direction,  Poisson's  ratio  ,  and 
the  anisotropy  ratios  need  to  be  defined.   Elasticity  matrix 
[D]  is  evaluated  using  equation  (4.41)  and  the  stiffness 
matrix   is  developed  following  the  procedure  described  in 
Appendix  A. 

For  nonlinear  properties  an  incremental  procedure 
is  required.   In  the  first  step  initial  values  for  Young's 
modulus  and  Poisson's  ratio  are  used.   In  subsequent  steps 
the  state  of  stress  due  to  all  previous  steps  are  read  from 
mass  storage  data.   To  determine  the  appropriate  values  of 
tangent  modulus  and  Poisson's  ratio,  the  interval  of 
octahedral  normal  stress  and  the  interval  of  the  ratio  of 
octahedral  shear  stress  to  octahedral  shear  stress  at  failure 
must  be  available. 
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For  octahedral  normal  stress,  the  program  conducts  a 
search  in  a  table  generated  in  OVERLAY  (1,0)  and  stores  the 
limits  of  that  particular  interval.   Then  the  octahedral 
shear  stress  at  failure  for  a  given  octahedral  normal  stress 
is  determined  directly  from  the  test  data,  in  the  case  of 
plane  strain  tests  or  from  the  following  expression  in  the 
case  of  plane  stress  tests: 


oct 


oct 


l+4> 


failure 


6N 


2(4>   -  <f  +  1)  - 


(N«fc  +  1) 


(5.2) 


ty  = 


al  +  a3 


1+Sin(j) 
1-Sin<j) 


where 


t     =  octahedral  shear  stress  at  failure 
oct 

a  =  octahedral  normal  stress  at  failure 

oct 

4)  =  friction  angle  of  soil 

°1  =  major  principal  stress 

a~  =  intermediate  principal  stress 

o.  =  minor  principal  stress. 

Derivation  of  expression  (5.2)  can  be  found  in 

Appendix  D. 


values  of 


In  equation  (5.2),  for  given  values  of  \\>   and  <j>  , 

T  _ 

can  be  determined.   Multiplying  this 


oct 


oct 


ratio  by  a   ..  the  value  of  t   .at  failure  is  known, 
1      oct  oct 
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the  failure  stress 
f 


Dividing  x  „.  of  the  element  by  t 

oct  J      oct 

ratio  can  be  determined.   This  ratio  is  the  second  parameter 
required  to  find  appropriate  value  of  tangent  modulus  and 
Poisson's  ratio;  the  program  automatically  interpolates  the 
stored  data  to  obtain  values  corresponding  to  the  actual 
state  of  stress. 

The  elasticity  matrix  [D]  for  the  soil  element  is 
calculated  using  equation  (4.41).   The  rest  of  the  proced- 
ure is  same  as  described  in  Appendix  A. 

Structural  Assembly 

Knowing  the  stiffness  matrix  of  an  element,  the 
structure  can  be  assembled  by  a  routine  operation  of  adding 
appropriate  coefficients  in  the  total  stiffness  matrix  using 
element  information  supplied  by  OVERLAY  (1,0).   In  such  a 
way  the  first  block  of  equations  are  formed.   As  the  stiff- 
ness matrix  is  symmetric  in  nature,  only  the  above  diagonal 
portion  of  the  matrix  is  generated  and  stored  in  a  random 
mass  storage  unit.   The  same  memory  space  can  now  be  used  to 
assemble  the  second  and  subsequent  blocks  of  equations  until 
all  nodes  are  considered. 

To  calculate  the  nodal  load  vector  due  to  gravity 
load,  total  weight  of  each  soil  element  is  determined  by 
multiplying  unit  weight  by  the  area.   Total  load  of  an  ele- 
ment is  distributed  in  all  six  nodes.   One-twelfth  of  total 
load  is  carried  by  each  of  corner  node  and  one-quarter  of 
total  load  is  carried  by  each  mid-side  node.   Total  load  at 
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a  node  is  the  summation  of  nodal  loads  from  each  element 
connected  at  the  node.   This  nodal  load  vector  is  also 
stored  in  a  random  mass  storage  unit. 

Solution  of  Stiffness  Equation: 

Several  approaches  for  solution  of  simultaneous 
linear  equations  have  been  tried,  including  Cholesky's 
method,  and  modified  Cholesky's  scheme.   These  schemes  take 
advantages  of  banded  nature  and  symmetry  of  equations.   The 
Cholesky's  scheme  is  difficult  to  use  in  large  problems  be- 
cause the  method  requires  the  complete  matrix  to  remain  in 
central  memory  during  the  solution  process,  and  requires 
taking  square  root.   The  square  root  of  a  negative  argument 
is  imaginary,  so  the  scheme  requires  a  very  well  conditioned 
matrix  (which  may  not  be  true  for  conditions  like  slip  in 
the  interface  element  or  very  low  tangent  modulus  value  at 
stress  ratio   1)  .   This  difficulty  can  be  overcome  by  using 
modified  Cholesky's  method  which  avoids  taking  square  root, 
but  for  large  problems  the  computer  time  needed  for  solution 
is  prohibitive. 

Christian  (1973)  reviewed  the  existing  schemes  for 
solution  of  simultaneous  equations.  He  pointed  out  advan- 
tages and  disadvantages  of  different  methods.  But  none  of  the 
methods  reviewed  can  be  considered  efficient  (and  economic) 
under  all  conditions.  Klyuyev  and  Kokovkin-Scherbak  (1965) 
pointed  out  that  for  solution  of  a  system  of  equations  the 
Gauss  elimination  process  requires  the  least  number  of 
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arithmetic  operations.   Based  on  this  finding  a  scheme  for 
solution  has  been  developed  using  Gauss  elimination.   In  this 
process,  the  equations  are  stored  on  disc  unit  sequencially . 
During  the  reduction  process  only  a  block  of  equations  of 
the  size  of  a  semi-band  width  by  semi-band  width  remains  in 
central  memory.   After  reduction,  one  block  of  reduced 
equations  is  written  on  a  tape  unit  and  one  block  of  equation 
is  read.   Solution  is  obtained  by  backsubstitution  in  the 
reduced  equation.   This  is  done  in  a  similar  fashion  as  in 
the  reduction  process  itself.   Solution  of  the  stiffness 
equations  for  given  nodal  load  and  boundary  condition  is  per- 
formed in  OVERLAY  (6,0)  whose  flow-chart  is  shown  in  Figure 
5.8. 

From  the  node  point  deflections,  the  stresses  and 
strains  are  determined.   This  is  done  in  OVERLAY  (4.0)  and 
is  shown  in  Figure  5.6.   Different  components  of  strain  can 
be  determined  from  geometric  properties,  shape  functions  and 
nodal  displacements  by  using  equation  (3.12).   Once  strains 
are  known,  the  stresses  can  be  determined  by  multiplying  the 
strain  matrix  with  [D]  matrix  for  the  same  element.   Stress 
and  strain  increments  for  each  increment  cycle  are  evaluated 
and  added  to  the  values  of  the  previous  cycle. 

The  last  and  final  step  in  the  scheme  is  to  print  the 
results.   This  is  done  in  OVERLAY  (5.0)  and  shown  in  Figure 
5.7.   In  this  OVERLAY  the  node  point  deflections  are  added 
to  those  from  previous  cycles  and  all  components  of  element 
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FIND  NODE  POINT  DISPLACEMENTS  FOR 
AN  ELEMENT 


COMPUTE  STRAINS  AND  ADD  TO  RESULTS 
OF  PREVIOUS  CYCLE 


COMPUTE  STRESSES  AND  ADD  TO 
RESULTS  OF  PREVIOUS  CYCLE 


ALL 
ELEMENTS 
CONSIDERED? 
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Figure  5.6.   Flow  Chart  for  OVERLAY  (4,0) 
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Figure  5.7  Flow  Chart  for  OVERLAY  (5,0) 
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READ  STIFFNESS  MATRIX  AND  LOAD  IN 
BLOCKS  FROM  MASS -STORAGE  UNIT 


MODIFY  STIFFNESS  MATRIX  FOR  BOUNDARY 
CONDITIONS 


REDUCE  STIFFNESS  MATRIX  AND  LOAD 
VECTOR,  STORE  ON  TAPE 


SOLVE  BY  BACKSUBSTITUTION  AND  FIND 
NODAL  DISPLACEMENTS 


Figure  5.8.   Flow  Chart  for  OVERLAY  (6,0) 
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stress  and  strains  and  displacements  are  printed. 

In  incremental  analysis  OVERLAYS  (4.0),  (6,0)  and 
(5,0)  are  repeated  for  the  desired  number  of  load  incre- 
ments.  For  analysis  of  sequence  of  construction,  in  steps, 
the  height  of  each  layer  has  to  be  defined.   According  to 
the  given  height,  the  program  will  select  which  elements 
will  appear  in  a  particular  layer.   Accordingly,  nodal  loads 
due  to  self-weight  for  those  elements  only  will  be  considered 
for  analysis.   This  follows  the  routine  procedure  of  assembly 
and  solutions  described  above. 

No-tension  Analysis 

It  is  well-known  that  most  soils  have  very  poor 
resistance  to  tensile  force,  and  in  most  cases  can  be 
assumed  to  have  no  resistance  in  tension.   For  this  purpose 
the  analysis  has  to  be  modified  so  that  tension  is  reduced 
to  zero  (or  to  a  specified  small  value).   In  other  words, 
the  tensile  forces  have  to  be  eliminated  from  those  elements 
which  experience  tension  and  the  effect  of  this  stress  re- 
lease in  the  adjoining  elements  has  to  be  evaluated.   Unfor- 
tunately, very  limited  research  has  been  done  in  this  partic- 
ular area.   Zienkiewicz,  et.  al.  (1968)  developed  a  concept 
for  'no-tension'  analysis.   The  major  steps  in  the  analysis 
are  summarized  below: 

(1)   Determine  those  elements  in  which  tension 
exists . 
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(2)  Calculate  the  nodal  loads  that  produced  the  ten- 
sile stresses. 

(3)  Apply  equal  and  opposite  nodal  forces  in  each 
element  where  tension  exists  and  recalculate  the  stresses  in 
all  the  elements. 

(4)  If,  at  the  end  of  step  (3) ,  tensile  stresses 
are  still  present  steps  (2)  and  (3)  are  repeated  until  con- 
vergence is  reached. 

Chang  and  Nair  (1973)  reported  that  use  of  this  tech- 
nique requires  a  large  number  of  iterations  and  in  some 
cases  the  solution  may  never  converge.   The  reason  for  this 
uncertainty  lies  in  the  fact  that  there  is  no  rational 
basis  for  evaluating  equivalent  nodal  loads  in  step  (3)  for 
a  two-dimensional  state  of  stress.   Chang  and  Nair  (1973) 
proposed  a  modification  by  altering  the  value  of  Poisson's 
ratio,  which  may  increase  the  rate  of  convergence  but  does 
not  alter  the  uncertainty  of  convergence. 

Development  of  tensile  stresses  in  interaction  ele- 
ment causes  a  similar  problem.   Also  when  the  ratio  of  shear 
stress  to  normal  stress  in  the  interaction  element  exceeds 
the  limiting  value,  excess  shear  stress  also  has  to  be  elim- 
inated.  A  method  for  estimating  nodal  loads  to  eliminate 
excess  tension  or  shear  is  proposed,  considering  the  follow- 
ing aspects : 

(1)   The  state  of  stress  before  application  of  the 
load . 
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(2)  Increments  of  stresses  due  to  the  current  load 
increment. 

(3)  In  this  manner  the  stiffness  matrix  is  kept  un- 
changed, which  eliminates  the  need  for  reevaluation  of  ele- 
ment stiffnesses  and  then  assembly  to  form  the  complete 
structure . 

The  procedure  for  no-tension  analysis  adopted  in  this 
study  is  described  below.   Figure   5.9  shows  a  state  of 
stress  o.  in  an  element,  which  represents  the  summation  of 
stresses  up  to  ith  increment  of  load.   Aa .  ,  represents 
change  in  stress  due  to  (i+1)  th  load  increment  and  a .  ,  is 
the  total  stress  at  the  end  of  this  increment  which  is  a 
net  tension  in  the  element.   The  tensile  stress,  a . ,,  has 
to  be  eliminated  by  applying  equivalent  nodal  loads  on  each 
node  of  the  element.   Equivalent  nodal  load  for  an  element 
is  given  by  the  following  expression 

{P}  -  [k]  •  {du}  (5.3) 

where 

{P}  =  Equivalent  nodal  load  vector 

[k]  =  Element  stiffness  matrix 
{du}  =  Nodal  displacement  vector. 

Equivalent  nodal  load  for  this  element  is  evaluated 
using  equation  (5.3)  where  {du}  is  the  nodal  displacement 
vector  for  (i+1)  th  increment  of  load.   The  load  increment 
is  modified  by  multiplying  by  a  load  factor  F,  which  is 
given  by  the  following  equation: 
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FIGURE    5.9    SCHEME    FOR   'NO -TENSION'    ANALYSIS 
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\oL\ 

F  =  1  - ,  where     denotes  absolute  value 

Aa 
1   i+l' 

(5.4) 

Then  the  correcting  nodal  load  is  given  by 

{P)c  =  -  F  •  {P}  (5.5) 

where 

{p }   =  Correcting  load  vector, 
c 

The  correcting  nodal  loads  (P)   are  applied  to  all 
nodes  of  an  element  which  developed  tension.   A  similar  pro- 
cedure is  applied  to  all  elements  in  tension.   The  complete 
structure  is  solved  again  with  only  the  correcting  loads 
with  the  structure  stiffness  matrix  kept  unchanged.   In- 
crements in  stress  are  added  to  those   from  previous  incre- 
ment and  checked  for  tension.   If  the  tension  is  within  a 
small  acceptable  limit,  the  next  load  increment  is  applied. 
On  the  other  hand,  if  tension  exists,  the  procedure  of 
stress  release  is  repeated  for  a  specified  number  of  times. 
If  the  procedure  fails  to  converge  after  the  specified 
number  of  iterations,  further  execution  is  stopped. 

Development  of  tension,  or  of  stress  ratio   in  ex- 
cess of  the  specified  limit,  in  the  interaction  element 
causes  similar  problems,  and  a  similar  procedure  is  applied 
in  such  cases.   The  load  factor,  F     in  equation  (5.  5)  is 
determined  in  the  following  manner. 
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Let  t  and  a  be  the  total  shear  and  normal  stresses 
in  an  interaction  element  at  the  end  of  ith  load  increment. 
At  and  Aa  are  increments  of  x  and  a  due  to  (i+l)th 
increment. 


x 
R-    =  —  =  stress  ratio  after  ith  load  increment. 

1     a 

t+ At 
Ri+1  =  a+Aa~  ~  stress  ratio  after  (i+1)  load  increment. 

R.  ,  is   greater  than  limiting  value  of  stress  ratio,  R^ . 

_L  i  X  J_j 

At 
Also    let  Q   =   t—  for    (i+l)th   increment. 

a • At    -   t    • Aa 

AR  =    RT     -    R,       =   -. — ~-  (5.6) 

h  1  0' (o   +   Aa    ) 


or 


Aa       (a«Q   -   t) 
AR     =   or-\o   +   Aar)  (5'7) 

Solving  for  Aa   from  (5.7) 

.  AR«a 

Aa   = 


r    (Q  -  R1  -  AR)  (5>8) 

where  Aa   is  the  increment  of  a  required  to  maintain  the 
stress  ratio  of  the  limiting  value. 

For  release  of  excess  shear,  i.e.,  to  maintain  stress 
ratio  within  allowable  limit,  equivalent  nodal  loads  given  by 
equation  (5.5)  has  to  be  applied  in  which  F  is  given  as 

Aar 

F  =  1.  -  T-i  (5.9) 

Aa 

andAo   is  defined  in  equation  (5.8). 
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This  procedure  also  suffers  from  uncertainty  of  con- 
vergence, and  further  research  is  needed  in  this  area. 
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CHAPTER  VI 
FINITE  ELEMENT  ANALYSIS 

6 . 1   Selection  of  Boundaries  and  Boundary  Conditions. 

To  analyze  the  culvert  problem  by  the  finite  element 
method,  a  set  of  finite  boundaries  must  be  established  and 
the  conditions  at  these  boundaries  defined  realistically. 
A  decision  must  be  taken  to  establish  (a)  the  size  of  the 
region  to  be  considered,  (b)  the  size  of  elements  in  differ- 
ent parts  of  the  region,  and  (c)  what  boundary  conditions 
should  be  used  at  each  boundary.   Depending  on  the  type  of 
problem,  type  of  analysis,  major  points  of  interest,  time 
and  resources  available  for  the  study,  the  answer  to  the 
above  questions  may  vary  over  a  wide  range.   However,  in- 
sight into  the  problem,  past  experience,  and  some  trial 
solutions  help  in  the  decision  making  process. 

In  the  usual  case,  symmetry  across  a  vertical  axis 
through  the  center  of  culvert  may  be  assumed  except  for 
special  cases  like  nonsymmetric  external  loading,  non- 
uniform backfilling  sloping  surfaces  and  similar  situations. 
The  assumption  of  symmetry  across  the  axis  reduces  the 
problem  to  one-half  size.   Trial  runs  revealed  that  the 
lateral  (vertical)  boundary  should  be  at  a  distance  three 
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to  four  times  the  radius  of  pipe,  and  the  bottom  (horizon- 
tal) boundary  at  least  one  diameter  below  the  pipe  invert. 

Selection  of  the  same  size  of  element  throughout  the 
medium  is  uneconomical  both  in  terms  of  time  required  to 
digitize  the  problem  and  cost  of  run.   The  general  practice 
is  to  use  smaller  sized  elements  in  zones  of  primary  inter- 
est and  in  areas  of  relatively  high  stress  gradients.   Large 
elements  can  be  used  elsewhere,  but  the  transition  between 
large  and  small  elements  should  be  gradual.   In  the  culvert 
problem,  the  major  area  of  interest  is  in  the  zone  near  the 
pipe,  hence  smaller-sized  elements  are  used  in  a  zone  of  half 
the  radius  of  pipe  and  the  size  increased  gradually  towards 
the  boundaries.   For  numerical  stability,  "flat"  triangles 
are  avoided. 

Conditions  at  boundaries  may  be  defined  either  by 
defining  nodal  forces  or  displacements.   Both  approaches  are 
acceptable  but  it  is  generally  easier  to  visualize  boundary 
conditions  in  terms  of  displacements  rather  than  nodal 
forces.   The  boundary  conditions  used  in  this  study  are  as 
given  below. 

(1)  Vertical  Axis  of  symmetry  -  Zero  horizontal  dis- 
placement; free  movement  vertically. 

(2)  Horizontal  Lower  Boundary  -  Zero  vertical  and 
zero  horizontal  displacement. 

(3)  Vertical  boundary  laterally  from  center  of  pipe 
-  Conditions  at  this  boundary  are  the  same  as  those  for  the 
axis  of  symmetry. 
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These  assumptions  of  boundary  conditions  are  reason- 
able and  are  commonly  used  (Corotis,  Farzin,  Krizek,  1974). 

6 . 2   Selection  of  Example  Problems 

The  main  purpose  of  this  study  is  to  develop  an 
analytical  tool  which  can  adequately  predict  the  performance 
of  culverts  buried  in  soil  under  various  situations.   The 
scheme  for  analysis  and  the  computer  code  developed  has  been 
checked  for  internal  consistancy  using  problems  for  which 
known  solutions  exist.   To  examine  the  capability  of  this 
program  to  analyze  culvert  problems  in  a  realistic  fashion 
a  pipe  culvert  of  varying  stiffness  was  chosen  and  its  be- 
havior analyzed  for  different  soil  conditions,  various  con- 
struction procedures  and  different  amounts  of  limiting  slip 
at  the  soil-pipe  interface.   A  corrugated  aluminum  pipe 
was  selected  as  an  example  problem  for  study.   The  pipe  has 
a  nominal  diameter  at  10  feet  and  is  embedded  in  soil  with 
20  feet  (or  40  feet)  height  of  cover  (Figure  6.1).   The  prob- 
lem has  been  idealized  using  340  nodes  and  174  elements 
(Figure  6.2).   The  lower  boundary  has  been  fixed  at  10  feet 
below  the  bottom  of  the  pipe  and  the  lateral  boundary  has 
been  chosen  at  a  distance  of  32  feet  from  the  center  line  of 
pipe. 

Properties  of  Soil 

The  backfill  material  selected  is  Monterey  No.  0 
sand  which  is  mainly  composed  of  quartz  and  feldspar  grains. 


Top  of    Embankment 


"liSF- 


Corrugated     Aluminum 
Pips 


A=  0.146   ir#in 

I  =  0.136  in4/in 
E  =    I07  psi 


H- 


=   0.33 


.37 


-T7JWT 


20'    OR    40' 


10' 


Sand    Backfill 

Unit    Wt.  -120  pcf  10 

E     =l.55xl05psf  I 

JUL    =  0.40 

i  i  >  )  i  i  r~rrr !  / 1  n  1 1  n  >  t~t  >  i  i  i    >  i    >  >  >   A)    >   i  >   /  -> — 7— t—t 


Hard    Clay 


FIGURE    G.I     GEOMETRY    OF   SAMPLE    CULVERT 
PROBLEM 
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FINITE  ELEMENT  MESH 


>- 


X  Ft. 


FIGURE  6.2    FINITE    ELEMENT    MESH     FOR    EXAMPLE 
PROBLEM. 
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It  is  a  uniform  medium  sand  with  uniformity  coefficient  C 

2  u 

=  d,A/d,_  =  1.53,  and  an  average  diameter  of  dCA  =  0.43  mm. 

Specific  gravity  is  2.645.   Minimum  and  maximum  void  ratios 

are  e  .   =  0.656,  e    =  0.860.   From  triaxial  tests,  the 
min  max 

friction  angle  <J>  for  the  sand  is  34.6  degrees  when  e  =  0.78. 
Lade  (1972)  conducted  plane  strain  tests  on  rectang- 
ular prismatic  specimens  10.6  cm  high,  4.6  cm  wide  and  11.5 

2 
cm  long  confining  pressures  a~   of  0.3,  0.6  and  1.2  kg/cm 

in  drained  conditions.   Stress-strain  and  volume  change 
characteristics  obtained  are  shown  in  Figure  4.8.   These 
results  have  been  used  to  represent  nonlinear  stress  depend- 
ent soil  properties  (Figures  4.9  and  4.10). 

Finite  Element  Idealization 

Sixteen  curved  bar  elements  were  utilized  to  idealize 
the  pipe,  i.e.,  the  angle  made  by  each  segment  of  pipe  at  the 
center  is  11.25  degrees.   The  total  number  of  interaction 
elements  is  16.   340  nodes  and  174  elements  have  been  used 
to  idealize  the  problem.   The  problem  being  symmetric,  only 
a  half-section  passing  through  center-line  of  pipe  was 
analyzed.   The  finite  element  mesh  is  shown  in  Figure  6.2. 

Problems  Selected 

In  order  to  study  effects  of  different  material 
parameters  and  cf  considering  construction  layers  with  and 
without  allowing  slip  to  occur  at  the  soil-pipe  interface,  a 
variety  of  problems  was  analyzed. 
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Pipe  Properties.  Corrugated  aluminum  nipe 

EI  =  1.36  x  10  lb.  in  /inch 

E  =  10T  lb/in2 

V  =   0.33 

I  =  0.136  in  /in 

Radius  R  =  60  inch 

2 
Area  A   =  0.1^6  in  /inch 

Soil  Properties.   Moderately  loose  sand 

Unit  weight  y     =  120  Ibs/cu.  ft. 

m 

(Elastic)  Tangent  Modulus,  E  =  1.55  x  105  psf. 
(Elastic)  Poisson's  Ratio,  jj  =  OJ4O 
The  effect  of  varying  E  was  examined 
Nonlinear  properties  =  direct  input  of  plane  strain 
test  data.   The  effect  of  soft  inclusions  shown  as 
a  speckled  zone  in  Figure  6.2  was  analyzed. 
The  construction  sequence  was  broadly  classified  into  the 
following  categories: 

I.   Single  Layer  -  Culvert  and  backfill  are  placed 
in  a  single  operation  with  the  weight  applied 
either  as  a  single  load,  or  in  increments. 
II.   Multi-layer  -  In  this  case  the  pipe  and  backfill 
are  constructed  in  layers.   The  layer  heirhts  and 
number  of  load  increments  tier  layer  are  shown  in 
Table  6.1.' 
The  height  of  backfill  was  also  varied,  with  20  feet 
and  hO   feet  heights  of   cover  considered.   Table  6.2  lists  the 
range  in  parameters  analysed  for  a  single  layer  of  backfill. 
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Table  6.1(a)   Layer  Heights  for  Incremental  Analysis,  20  Ft. 
Fill. 


Layer 
No. 

Layer  Height   (Ft. ) 

No.  of  Increments 

From 

To 

A 

B 

1 

0 

9.5 

1 

2 

2 

9.5 

12.0 

1 

2 

3 

12.0 

15.0 

1 

2 

4 

15.0 

18.0 

1 

2 

5 

18.0 

21.0 

1 

2 

6 

21.0 

26.0 

1 

2 

7 

26.0 

32.0 

2 

3 

8 

32.0 

40.0 

2 

3 

Total 

10 

18 

Table  6.1(b)   Layer  Heights  for  Incremental  Analysis,  40  Ft. 
Fill. 


Layer 
No. 

Layer  Height   (Ft. ) 

No.  of  Increments 

From 

To 

1 

0 

9.5 

2 

2 

9.5 

12.0 

2 

3 

12.0 

15.0 

2 

4 

15.0 

18.0 

2 

5 

18.0 

22.0 

2 

6 

22.0 

32.0 

2 

7 

32.0 

44.0 

6 

8 

44.0 

60.0 

6 

Total 

24 

1^2 


Table  6.2 
Example  Problem  Description  -  I.  Single  Layer,  20  Feet  of  Cover 


No-        ESoil  EIPiue  No*  of      Interface 

Increments     Slin 


L-l  E±  =   1.55xl05  psf     EI  =   1.13xl05   lb   ft2/ft  1  No 


v.  =  O.k 

l 

V  =  0.33 

L-2 

*Ji 

EI/2 

L-3 

Ei 

EI/5 

L-L 

E. 

l 

EI/10 

L-5 

E. 

l 

2EI 

L-6 

Ei 

5EI 

L-T 

E. 
l 

10EI 

L-8 

E. 

l 

20EI 

L-9 

E. 

l 

50EI 

L-10 

i 

100EI 

L-11 

E./10 

EI/10 

L-12 

E1/10 

EI/5 

L-13 

E./10 

EI/2 

L-ll» 

E./10 

EI 

L-l  5 

Ej/10 

2EI 

L-l  6 

E./10 

5EI 

L-17 

E./10 

l 

10EI 

L-18 

E./10 
i 

20EI 

L-19 

E./10 

50EI 

L-20 

E.  /10 

100EI 

L-21 

10E. 

l 

EI/10 

1 

No 

1 

No 

1 

No 

1 

No 

1 

No 

1 

No 

1 

No 

1 

No 

1 

No 

1 

No 

1 

No 

1 

No 

1 

No 

1 

No 

1 

No 

1 

No 

1 

No 

1 

No 

1 

No 

1 

No 

Table  6.2  (Continued) 
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No. 


JSoil 


EI 


L-22 

10E. 

1 

L-23 

10E. 

l 

L-2U 

10E 

L-25 

10E 

L-26 

10E. 

l 

L-27 

10E. 

l 

L-28 

10E 

L-29 

10E. 

l 

L-30 

10E. 

l 

L-31 

Ei 

L-32 

Ei 

L-33 

Ei 

L-31** 

E. 

E, 


E 


Pipe 


No.  of 
Increments 


Interface 
Slip 


L  -  Linear  elastic  soil  properties 


EI/5 

1 

No 

EI/2 

1 

No 

EI 

1 

No 

2EI 

1 

No 

5EI 

1 

No 

10EI 

1 

No 

20EI 

1 

No 

50EI 

1 

No 

100EI 

1 

No 

EI/100 

1 

No 

EI/50 

1 

No 

EI/20 

1 

No 

EI 

1 

No 

EI/10 

1 

No 

JMEI 

1 

No 

lUEI 

1 

No 

100EI 

1 

No 

10EI 

1 

No 

E  =  1.0x10  psf,  v  =  0.47 
J  j 


E„„,-i  ~  77~  near  springline 
soil    ioo 


1*0  ft.  height  of  cover 
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Table  6.2  (Continued) 


No. 

ESoil 

EIPipe 

No.  of 
Increments 

Interface 
Slip 

NL-1 

NL 

iHei 

5 

No 

NL-2 

NL 

lJ+l+EI 

5 

Yes 

NL-3 

NL 

ikkEl 

10 

No 

NL-1* 

NL 

ikkEl 

10 

Yes 

**NL-5 

NL 

iUei 

5 

No 

**NL-6 

NL 

iMei 

5 

Yes 

**NL-7 

NL 

lkMi 

5 

No 

**NL-8* 

NL 

IkkEl 

5 

No 

**NL-9* 

NL 

iUei 

5 

Yes 

Parameters  analyzed  using  8  layers  to  simulate  actual  construction  are 
in  Table  6.3.  The  number  of  increments  per  layer  are  as  shown  in  Table 
6.1(a)  and  6.1(b)  for  20  and  1+0  feet  of  cover,  respectively. 


Table  6.3 
Example  Problem  Description  -  II.  Multi-Layer,  Nonlinear  Soil,  20  Feet  of  Cover 


N°-         ESoil  EIP1r)„  Total  No.  of  Interface 


Pipe 


Increments  Slip 


ML-1         NL  1UEI  10  No 

^-2         NL  1UUEI         10  Yes 

ML-3         NL  il,l,EI 


18 


ML"1*         NL  U^EI 


No 


NL  IkkEl  l8  Yc: 


**ML-5 


NL  IkkEl  2k 


NL  -  Nonlinear  soil  properties,  d'irect  "in^ut  "of  TesTditV 


No 


*E 


soil  "  10  near  sPriKS  line;  »«1*0  feet  of  cover 
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CHAPTER  VII 
PRESENTATION  AND  DISCUSSION  OF  RESULTS 

The  primary  purpose  of  the  example  problems  is  to 
study  changes  in  response  of  the  pipe  to  different  input 
parameters,  and  to  recognize  which  of  the  variables  have 
significant  effects  on  the  behavior  of  pipe.   The  factors 
chosen  for  comparison  to  reflect  culvert  response  were: 

(1)  Deformed  shape  of  pipe  and  percent  change  in 
vertical  diameter,  A  Y% . 

(2)  Maximum  thrust  in  pipe,  P 

c   r  '   max 

(3)  Maximum  moment  in  pipe,  M 

r  r  '   max 

(4)  Maximum  extreme  fiber  stress  in  pipe,  ofl 

max 

(5)  Marston-Spangler  E'  value  back-calculated  from 
equation  (2.4) . 

(6)  Horizontal  and  vertical  stress  in  the  soil 
adjacent  to  the  springline  of  the  pipe,  a   and  a  ,  respec- 
tively. 

(7)  Vertical  stress  in  the  soil  at  the  crown  of  the 


pipe,  an. 


(8)   Maximum  value  of  ratio  of  shear  stress  to  normal 


stress  at  the  soil-pipe  interface,  — 


max 
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When  the  effects  of  slip  at  the  pipe-soil  interface  was  being  examined, 
the  limiting  value  of  -  was  set  at  0.33.   To  increase  the  rate  of  convergence 
a  10%  tolerance  in  this  ratio  was  allowed;  also,  a  small  amount  of  tensile 
stress  [not  exceeding  120  psf  ]  was  allowed  in  the  no-tension  analysis. 

Solutions  to  example  problems  were,  organized  into  three  groups: 

(1)  Problems  using  linear  elastic  soil  properties,  varying  the 
relative  stiffness  of  pipe  and  soil. 

(2)  Problems  using  non-linear  soil  properties  in  which  the  effects 
of  construction  sequence  are  analyzed. 

(3)  Problems  similar  to  (2)  but  using  a  larger  fill  height. 
Group  One : 

A  listing  of  problems  solved  in  the  first  group  is  given  in  Table  6.2. 
Problems  L-l  to  L-10  show  the  effect  of  changing  pipe  stiffness  keeping  the 
surrounding  soil  unchanged.   Problems  L-ll  to  L-12  and  L-21  to  L-31  cover  the 
same  range  in  pipe  stiffness  but  the  soil  modulus  is  reduced  by  a  factor  of 
10,  and  increased  by  a  factor  of  10,  respectively.   In  all  cases  the  problems 
were  solved  with  the  full  gravitational  load  due  to  self-weight  of  soil  applied 
in  one  operation.   A  summary  of  the  results  obtained  is  presented  in  Table  7.1. 
Percent  change  in  vertical  diameter,  AY%,  for  various  values  of  EI  are  plotted 
on  a  semilog  scale  and  shown  in  Figure  7.1a  for  each  soil  stiffness.   It  is 
seen  that  compression  of  the  pipe  in  the  vertical  direction  Is  dependent  on 
the  stiffness  of  pipe  and  on  the  modulus  value  of  the  surrounding  soil.   The 
soil  modulus  has  an  enormous  influence  on  the  deflection  of  flexible  pipes  but 
the  effect  becomes  less  and  less  pronounced  as  the  pipe  stiffness  increases. 
If  the  soil  and  pipe  properties  are  homogeneous  and  linear,  the  change  in 
pipe  diameter  is  significant  only  if  the  soil  is  weak  and  the  pipe  is  flex- 
ible. 
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Maximum  moment  and  maximum  thrust  in  the  pipe  are  plotted  vs.  pipe 
stiffness  in  Figure  7.1b  for  each  soil  stiffness.   It  is  seen  that  the 
maximum  moment  increases  sharply  with  increasing  pipe  stiffness — especially 
in  the  case  of  a  soft  soil- -but  the  maximum  thrust  in  the  pipe  is  virtually 
unaffected  [it  depends  primarily  on  height  of  fill  and  diameter  of  pipe]. 

The  maximum  extreme  fibre  stress  in  the  pipe  is  shown  plotted  vs.  pipe 
stiffness  in  Figure  7.1(c) .   Although  soil  stiffness  is  a  factor — with 
weaker  soils  inducing  higher  stresses — the  predominant  effect  is  due  to 
pipe  stiffness.   The  relative  importance  of  thrust  vs.  moment  on  the  maximum 
extreme  fibre  stress  is  shown  in  Figure  7.1(d).   For  a  typical  corrugated 
metal  pipe  (normalized  EI  =  l)  85  percent  of  the  maximum  stress  is  due  to 
moment  (vs.  15  percent  due  to  thrust)  in  a  weak  soil;  even  in  a  moderately 
good  soil  (E.  =l)  the  ratio  is  still  57  percent.   Thus,  very  small  moments 
in  flexible  pipes  can  cause  yielding  of  the  pipe  wall.   If  local  buckling 
does  not  develop,  the  soil  pressures  on  the  pipe  will  redistribute  to  re- 
duce the  induced  moments;  if  not,  there  is  danger  of  collapse.   Thus,  local 
buckling  may  be  viewed  as  a  key  performance  criterion  for  flexible  culverts. 
In  stiff  pipes,  the  induced  moments  dominate  pipe  performance. 

In  all  the  above  cases  the  Marston-Spangler  values  of  soil  modulus  E' 

were  back-calculated;  they  are  listed  in  Table  7-1  and  shown  plotted  in 

Figure  7.1e.   If  E'  were  a  soil  property  its  value  would  be  constant  for  a 

given  E   .,  and  a  given  pipe  diameter.   It  can  be  seen  that  for  weak  soils 
soil 

the  E'  values  are  strongly  dependent  on  pipe  stiffness  for  the  same  E   .-, 
[E1  is  also  greatly  dependent  on  the  height  of  cover].   Accordingly,  E' 
is  not  a  soil  property,  and  there  is  no  hope  of  empirically 
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NORMALIZED     PIPE    STIFFNESS,   El  pipe/l.l33xl05  lbft2/ft 

FIGURE  7-Kb)     EFFECT    OF    RELATIVE    STIFFNESS    OF    PIPE 

AND    SOIL    ON    MAXIMUM    MOMENT  (Mmax )  AND 
MAXIMUM    THRUST  (Pmax)    IN    THE    PIPE. 
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correlating  values  of  E'  that  vould  be  valid  for  a  vide  range  of  culvert 
problems • 

To  investigate  the  effect  of  poor  backpacking  near  the  springline 
two  problems  (L-3^,  L-35)  were  solved  where  very  loose  material  is  placed 
in  this  location  (see  Figure  6.2).   Results  are  given  in  Table  7-1  and 
compared  below. 

£Y        aQ  (max)     a     at  spring 
EI        %_  psi  x  103      psf  x  103 

E   .,  =  1;  uniform  1         0.82        15-3  2.25 

soil 

Esoil  =  ^  1/10°  f         1        U.03       59.8         0.36 
springline 

In  this  case  the  same  pipes  are  subjected  to  much  higher  extreme  fibre 

stresses  and  the  deflections  at  the  crown  are  increased  by  a  factor  of  5. 

Lack  of  confinement  of  the  pipe  near  the  springline  adversely  affects 

culvert  performance,  a  fact  which  is  well  known  from  experience. 

In  problem  L-35,  the  pipe  is  ten  times  more  flexible  than  in  problem 
L-3^,  but  the  reductions  in  vertical  diameter  are  approximately  the  same. 
This  is  explained  by  the  formation  of  a  soil  arch  over  the  pipe  which 
transfers  load  away  from  the  crown.   Formation  of  the  "arch"  is  also  re- 
flected in  an  increase  in  vertical  stress  in  the  soil  away  from  the 
springline.   The  formation  of  a  soil  arch  is  the  principal  reason  why  it 
is  possible  to  construct  'super-span'  corregated  metal  culverts  (Fisher, 
I969).   It  is  seen  that  considerable  insight  into  culvert  behavior  [but 
not  the  prediction  of  culvert  performance]  can  be  gained  from  simple 
elastic  analysis. 

Figures  7-2  through  7-^  show  the  distribution  of  parameters  reflect- 
ing pipe  response  around  the  circumference 
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ORIGINAL    SHAPE 


MOMENT 
I03  ft  lb/ft 


DEFORMED    SHAPE 


l 1 

0        O.ift 


FIGURE  7.2    DISTRIBUTION    OF   MOMENT    AND   DEFORMED 
SHAPE    OF    PIPE,  PROBLEM    L-14. 
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FIGURE    7.3     DISTRIBUTION    OF    NORMAL    AND    SHEAR 
FORCES    IN    PIPE,  PROBLEM    L-14. 
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SHEAR     STRESS 
I02  psf 


NORMAL   STRESS 


FIGURE  7.4    NORMAL    AND    SHEAR    STRESS    IN 

INTERACTION    LAYER,  PROBLEM  L-14. 
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of  the  pipe  for  problem  L-14.   In  Figure  7.2,  deformed  shape 
of  pipe  and  distribution  of  nodal  moments  in  the  pipe  are 
shown.   Distribution  of  thrust  and  shear  force  in  the  pipe 
are  shown  in  Figure  7.3.   Distribution  of  normal  and  shear 
stresses  in  interaction  layer  (soil-pipe  interface)  is 
plotted  in  Figure  7.4.   The  results  are  typical  for  problems 
in  Group  1,  except  for  L-34  through  L-39. 

Group  Two 

The  main  purpose  of  the  problems  in  Group  2  is  to 
study  the  effects  of  nonlinear  soil  properties  and  simulation 
of  construction  sequence.   Results  of  problems  NL-1  to  NL-4 

and  ML-1  to  ML-4  are  listed  in  Table  7.2.   Also  for  compari- 

the 
son,  results  from  problem  L- 36  (which  is  one  layer,  linear 

elastic  solution)  have  been  included  in  the  same  table,  be- 
cause Young's  modulus  and  Poisson's  ratio  for  the  soil  is  the 
same  as  the  initial  values  in  the  non-linear  analysis. 

In  problems  NL-1  to  NL-4  construction  is  completed 
in  a  single  layer  but  the  total  load  is  applied  in  5  (for 
NL-1  and  NL-2 >  and  10  (for  NL-3  and  NL-4)  increments.   (See 
Table  6.2).   Results  obtained  for  10  increments  of  load  are 
very  similar  to  that  obtained  with  5  increments,  which  shows 
that  convergence  was  reached  after  5  increments  of  load  were 
applied  to  a  single  soil  layer.   For  other  geometries  and 
degrees  of  non-linearity  convergence  may  not  be  so  rapid. 

Comparison  of  results  with  and  without  interface 
slip  shows  very  little  or  no  variation.   This  is  because 
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very  few  interaction  elements  have  reached  failure  conditions  correspond- 
ing to  slip,  and  these  by  only  a  small  margin. 

Comparison  of  results  between  linear  (L-3&)  and  non-linear  soil 
properties  (NL-1)  shows  relatively  large  differences  in  deflected  shape 
and  in  the  extreme  fibre  stresses.  Also,  the  results  from  single  layer 
(5  to  10  increment)  nonlinear  analysis  are  comparable  to  the  results 
from  multi-layer  (8  layer)  nonlinear  analysis  which  shows  that  for  a 
corrugated  metal  culvert,  10  feet  in  diameter,  with  20  foot  of  cover  con- 
sisting of  a  granular  soil  the  effects  of  considering  nonlinear  soil 
properties  are  very  significant  but  consideration  of  the  construction 
sequence  is  unimportant  as  long  as  it  is  symmetrical. 

Problems  ILL— 1  to  ML-4  have  8  construction  layers  (see  Table  6.3)  and 
were  solved  using  nonlinear  soil  properties  and  incremental  analysis. 
In  this  case  also  interaction  has  very  little  effect  on  the  behavior  of 
pipe.   The  total  number  of  load  increments  (Table  6.3)  has  an  influence 
on  the  results  but  the  differences  between  10  and  18  increments  are  not 
large. 

In  the  construction  process,  the  pipe  is  placed  on  a  prepared  bed- 
ding and  the  placement  of  backfill  material  and  compaction  proceeds  in 
layers.   It  is  observed  during  this  backpacking  operation  that  the  pipe  elongates 
in  the  vertical  direction,  and  that  any  untied  vertical  struts  will  topple. 
In  the  finite  element  simulation  of  construction  in  layers  the  same  phe- 
nomenon is  observed;  in  all  solutions 
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(ML-1  to  ML-4)  the  pipe  actually  elongated  during  the  placement  of 
4th  and  5th  layer  of  backfull  material. 

Figures  7.5  to  7.7  show  the  distribution  of  parameters  reflecting 
pipe  response  around  the  circumference  of  the  pipe  for  problem  NL-1. 
The  deformed  pipe  shape  and  distribution  of  moments  are  shown  in  Figure 
7.5.  Distribution  of  thrust  and  shear  force  in  the  pipe  are  plotted  in 
Figure  7.6.   Distribution  of  normal  and  shear  stress  in  interaction  layer 
(at  pipe  interface)  is  shown  in  Figure  7.7.   Similar  plots  for  multi- 
layer solutions  are  given  in  Figures  7.8  to  7.10. 

Group  Three 

Problems  in  group  three  were  selected  to  show  the  influence  of 
increased  height  of  fill.   In  this  case  the  height  of  cover  has  been  increased 
from  20  feet  to  40  feet.   Results  for  linear  (L-37)  ,  nonlinear  (NL-5  to 
NL-9)  and  multi-layer  (ML-5)  analyses  are  tabulated  in  Table  7.3.   The 
differences  between  elastic  and  nonlinear  solutions  are  very  prominent 
in  this  case,  but  the  influence  of  slip  at  the  soil-pipe  interface  is  of 
no  importance.  Decreasing  pipe  stiffness  [NL-6  vs.  NL-7]  causes  increased 
shortening  in  vertical  diameter,  lower  moment,  and  higher  extreme  fibre 
stresses.  Making  the  pipe  more  flexible  creates  an  arching  effect  which 
is  evidenced  by  reduced  vertical  stress  in  the  soil  near  crown  of  pipe 
and  increased  vertical  stress  at  the  springline. 
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Moment     at    Pipe 
IC?  ft.  lb. /ft. 


Deflected    Shape  of    Pjpe 
10  ft 


0  0.5 


FIGURE    7.5     DEFLECTED    SHAPE    AND    DISTRIBUTION 
OF   MOMENT    IN    PIPE,  NONLINEAR 
ANALYSIS.  PROBLEM    NL-I. 
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Shear    Force    in    Pipe 
10  lbs/ft 


Normal    Force   in    Pipe 
10    lbs/ft 


FIGURE    76    DISTRIBUTION    OF   NORMAL    AND   SHEAR 
FORCE   IN   PIPE,  NONLINEAR    ANALYSIS, 
PROBLEM     NL-I. 
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Shear    Stress   in   Interaction 
10    psf. 


Normal    Stress    in   Interaction 
IC?  psf 


FIGURE   7.7    DISTRIBUTION    OF    NORMAL    AND     SHEAR 
STRESS   IN   INTERACTION    LAYER,  NON- 
LINER    ANALYSIS,  PROBLEM    NL- I . 
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Moment   in    Pipe 

I03  ft.  lb/ft. 

A     ML-  l 
o     ML -3 


Deflected     Pipe   Shape 

, ,       io2ft 

O    0.5 

O      ML-  I 

•      ML- 3 


FIGURE    7.8     DEFLECTED    SHAPE    AND    DISTRIBUTION 
OF    MOMENT   IN    PIPE,  MULTL  LAYER 
ANALYSIS,  PROBLEM    ML-  I ,  ML-  3. 
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Shear    Force  in    Pipe 
4 

10   lbs/ft 

D     ML-  I,  ML-  3 


Normal    Force  in    Pipe 
10   lbs/ft 

A      ML- I, 
O     ML- 3 


FIGURE    79   DISTRIBUTION    OF   NORMAL    AND    SHEAR 
FORCE   IN   PIPE  .  MULTI   LAYER    ANALYSIS, 
PROBLEM     ML- I,  ML-3. 
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Shear    Stress    in    Interaction 
10    psf. 


A    ml-  l 
O     ML -3 


Normal    Stress  in    Interaction 
I03psf. 


A     ML- I,  ML-3 


FIGURE  7.10    DISTRIBUTION    OF   NORMAL    AND    SHEAR 
STRESS    IN   INTERACTION    LAYER  .MULTI- 
LAYER   ANALYSIS,  PROBLEM   ML -I,  ML- 3 
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The  problem  with  40  feet  fill  has  been  solved  using 
less  stiff  material  near  springline  of  pipe  (problems  NL-8, 
NL-9) .  This  caused  increased  distortion  of  pipe,  and  a  35% 
increase  in  the  extreme  fiber  stresses.  Again,  slip  at  the 
soil-pipe  interface  with  a  limiting  x/a  =  0.33  had  no  in- 
fluence on  the  results — which  means  little  slip  occurred 
between  soil  and  pipe.   But  in  all  cases  the  maximum  thrust 

on  pipe  remains  virtually  unaffected  by  different  modes  of 
analysis. 

The  problem  with  increased  height  of  fill  has  been 
analyzed  using  8  construction  layers  and  24  load  increments 
(ML-5) .   Results  are  shown  in  Table  7.3.   Similarity  of  re- 
sults between  nonlinear  analysis  using  single  layer  (NL-5) 
and  multi-layer  construction  (ML-5)  is  striking.   For  a 
culvert  beneath  high  fills,  the  mode  of  construction  may  not 
be  of  importance,  as  the  total  height  of  fill  has  a  dominant 
effect.   This  may  not  be  true  with  nonhomogeneous  or  aniso- 
tropic backfill  material. 
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CHAPTER  VIII 
SUMMARY 

8. 1   Results 

Three  significant  improvements  were  made  in  the 
models  used  to  represent  components  of  a  culvert  problem: 

1.  An  element  consisting  of  a  curved  bar  with  bend- 
ing resistance  at  the  nodes  was  introduced  to  represent  seg- 
ments of  the  pipe.   It  was  shown  that  for  flexible  pipes 

(flexibility  number  <  10  )  large  errors  are  introduced  if  con- 
ventional (triangular  or  rectangular)  elements  are  used — un- 
less the  pipe  is  divided  into  a  very  large  number  of  elements. 
For  the  purpose  of  this  study,  the  properties  of  the  pipe 
material  were  restricted  to  linear  elastic  behavior. 

2.  An  isoparametric  triangular  element  with  one 
curved  side  (to  fit  the  pipe) ,  linear  strain  distribution, 
and  three  intermediate  nodes  was  used  to  represent  the  soil. 
After  many  trials,  considering  a  variety  of  ways  to  repre- 
sent soil  properties,  it  was  decided  to  use  the  results  of 
actual  test  data  expressed  in  terms  of  octahedral  stress  and 
strain.   The  results  are  stored  to  memory  and  the  computer 
automatically  interpolates  the  appropriate  soil  properties 
for  any  element  at  a  particular  stress  level  and  stress  in- 
crement. 
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3.   A  new  element  was  introduced — called  the  'inter- 
action' element—to  represent  behavior  at  the  soil  pipe 
interface.   The  element  has  zero  thickness  (essentially  zero 
displacement  normal  to  the  pipe)  but  can  induce  slip  when  the 
ratio  of  shear  stress  to  normal  stress  exceeds  a  prescribed 
limit  (angle  of  wall  friction) . 

A  computer  program  was  written  incorporating  the 
features  described  above  which  can  handle  situations  in  which 
the  soil  and/or  interaction  elements  are  incapable  of  resist- 
ing tension  and  in  which  the  loads  due  to  self-weight  of  the 
soil  can  be  applied  in  convenient  layers  (to  simulate  the 
construction  sequence) .   Anisotropic  stress-strain  properties 
in  the  pipe  material  and  in  the  soil,  and  inclusions  or 
zones  of  material  with  different  properties,  was  accounted 
for. 

8. 2   Conclusions 

Some  example  problems  were  solved  to  demonstrate  the 
versatility  of  the  program  and  to  investigate  the  influence 
of  such  factors  as  non-linear  soil  properties,  relative 
stiffness  of  pipe  and  soil,  inclusion  of  weak  materials  near 
the  spring  line,  and  construction  procedures.   From  these  re- 
sults the  following  preliminary  conclusions  are  drawn: 

A.  In  the  case  of  conventional  circular,  corrugated 
metal  pipes  buried  in  granular  soil,  if  the  fill 
is  placed  symmetrically 
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a)  considerable  insight  into  culvert  behavior 
can  be  gained  even  if  the  soil  is  represented 
by  a  linear  elastic  material 

b)  the  maximum  circumferential  thrust  in  the 
pipe  depends  mainly  on  pipe  diameter  and 
height  of  fill;  the  stiffness  of  the  soil  or 
the  relative  stiffness  of  the  pipe  is  of 
little  influence 

c)  although  the  maximum  moment  in  the  pipe  de- 
creases as  the  pipe  stiffness  decreases  the 
extreme  fibre  stresses  increase  sharply. 

This  suggests  that  a  key  factor  in  the  per- 
formance of  corrugated  metal  pipes  is  ade- 
quate resistance  to  local  buckling;  in  the 
absence  of  buckling,  yielding  of  the  pipe 
walls  will  redistribute  the  soil  pressures 
to  insure  stability.   Slip  in  bolted  seams 
(which  reduces  the  thrust  in  the  pipe  with- 
out impairing  the  bending  resistance) ,   is 
also  beneficial. 
B.   The  use  of  non-linear  soil  prpoerties  is  essen- 
tial to  the  prediction  of  culvert  performance, 
while  simulating  construction  sequence  (if  sym- 
metrical) is  of  lesser  importance — especially  as 
the  height  of  fill  increases  (to  40  feet  or  more) 

C.   Allowance  for  reduced  friction  at  the  pipe-soil  in- 
terface is  of  little  importance  for  10  foot  diameter 
corrugated  metal  pipes  buried  in  granular  fill 
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with  20  to  40  feet  of  soil  cover.   It  is  believed, 
however,  that  interface  slip  can  play  an  impor- 
tant role  in  the  behavior  of  large-span  corru- 
gated metal  arches. 
D.   The  Marston-Spangler  soil  modulus  E1  is  not  a 

soil  property — its  value  was  back-calculated  in 
the  example  problems  solved  and  was  shown  to  vary 
strongly  with  pipe  stiffness,  and 
amount  of  soil  cover.   Accordingly,  empirical 
correlations  of  E'  with  pipe  performance  are 
tenuous  and  cannot  be  used  as  a  rational  basis 
for  future  predictions. 

8.3   Recommendations  for  Future  Research 


Based  on  the  present  study,  it  is  recommended  that 
further  research  should  be  directed  in  the  following  areas. 

(1)   More  study  is  required  in  the  area  of  'no- 
tension'  analysis  and  in  the  allowances  made  for  interface 
slip.   A  methodology  is  needed  to  reduce  the  number  of 
iterations  required  to  reach  convergence. 

(2)   The  displacement  function  for  interaction  is 
linear,  for  soil  it  is  quadratic  while  that  for  pipe  is 
cubic.   Accordingly,  there  is  incompatibility  of  displace- 
ments along  the  element  interfaces  between  the  nodes.   The 
consequences  of  this  incompatibility  may  be  important  and 
should  be  investigated  further. 
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(3)  The  solution  scheme  needs  additional  study  to 
see  if  further  economies  can  be  effected  in  the  solution  of 
large  and  complex  problems,  possibly  including  automatic  gen- 
eration of  the  finite  element  mesh. 

(4)  The  development  of  plastic  hinges  and  of  local 
yielding  of  the  pipe  should  be  accounted  for.   This  is  neces- 
sary in  realistic  predictions  of  culvert  performance,  because 
if  yielding  is  allowed  (cracking  in  the  case  of  concrete 
pipes) ,  it  plays  an  important  role  in  the  economy  of  the 
structure. 

(5)  Sensitivity  studies  should  be  made  for  other 
pipe  geometries,  especially  for  the  pipe-arch. 

(6)  To  investigate  cases  in  which  the  differential 
settlements  in  the  longitudinal  direction  are  significant, 
a  three-dimensional  finite  element  analysis  should  be 
developed. 
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APPENDIX  -  A 
ISOPARAMETRIC  LINEAR-STRAIN  TRIANGULAR 
ELEMENT   WITH  ONE  CURVED  BOUNDARY 

Derivation  of  isoparametric  linear  strain  triangular 
finite  elements  has  been  presented  in  section  3.5  in  brief. 

A  more  detailed  treatment  is  described  here.   Figure 
A'l  shows  an  element  with  corner  nodes  marked  1,  2,  3  and 
intermediate  midpoint  nodes  are  4,  5,  6.   P(x,y)  is  any 
floating  point  inside  the  triangular  area.   Side  1-4-2  may 
be  curved.   The  triangle  may  be  subdivided  into  four  areas 

P23   =   area  A, 

P31   =   area  A_ 

1P2  and  the  curved  side   =   area  A~ 

Side  12  and  curved  side   =   area  A. 

Total  area  of  triangle  123   =   A* 

A*   =   A,  +  A2  +  A3  +  A4  (A.l) 

The  area  ratios  are 

?1      =      A*    '       ^2      =      A*         ^3      =      A* 
and  h   +    ^2    +    h      =      X      '      2t1  (A-2) 


where 


A4 
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Y  ♦ 
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□      Corner    Nodes 

V     Intermediate    Mid   Point     Nodes 


■*r 


FIGURE   A- 1       ISOPARAMETRIC    LINEAR    STRAIN 
TRIANGULAR    ELEMENT. 
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Relations  between   C-   and   x,y   system  are  given  in 
equation  (3.2).   Displacement  functions  F..  and  F.  have  been 
derived  in  equation  (3.5)  and  (3.6)  respectively.   Expres- 
sions for  F»,  F^,  F,.  and  F,  follows  the  same  procedure  and 
are  summarized  in  equation  (3.7)  and  shown  here. 


Fl   "  h 


=     K. 


I^l"1 


l-4n 


=  5- 


2?3-  i+4n" 


,.  ,  .2  €152 

(l-2n) 

-A p  r 

(l-4n)  ^2^3 
4 


(l-4n)  ^3^1 


i-4n 


(A. 3) 


In  vector  notation 

rn 

f   =  (f  f  f  f  f  f) 

-         1   2   3   4   5   6 

where  T  refers  to  transpose  of  matrix. 

Displacement  at  P  (Figure  A*l)  can  be  written  in 
terms  of  nodal  displacement  in  the  following  manner 


f   = 


F   0_ 
0   F 


u 

V 


(A. 4) 
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where   u   and   v   are  nodal  displacement  vectors  in   x  and 
y  directions,  respectively. 

Strain  in  the  element  can  be  determined  by  taking 
the  partial  derivative  of  equation  (A. 3)  with  respect  to   x, 
and   y.   Chain  rule  for  partial  derivatives  will  be  used, 
which  is 

5f      3   6F      5£. 


-  I 


1 


or 


6x      L      6£.       6x 
i=l   1 


5f      5   5F.    6£,    6F.   ££.,    5F.   6£_ 

_=  =  y  — 1  .  — L  +  — i  — 2   — 1  — 3 
i=l   x 


From  (A. 3) 
6F 


5? 

1 


~     =     I=2n  {4^i  "  X  +  2n} 


6F2       x 


5£2      l-4n  L  ""2 


6F3 


{4Co  "  1} 


6?3 

l^t>3  " 

5F4 

4C2 

eCi 

(i-2n)2 

5F4 

4^l 

652 

(l-2n)2 

(A. 5) 
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^5     4__ 

6?3     (l-4n)C2 


^5       4 p 

6£2      (l-4n)^3 


6F6 


6?!      (l-4n)s3 


6?3      (l-4n)^l 


65,        65. 


To  determine  derivatives  like   7 — -  and  -r— —  , 

Sx        Sy 

equation  (3.2a)  will  be  used  and  the  derivatives  can  be 


written  as 

6h 

= 

1 

6x 

2 (Area) 

Y2 

6C2 

= 

1 

6x 

2 (Area) 

"  Y31 

6C3 

= 

1 

6x 

2 (Area) 

Yl 

where   y . .   =   y .  -  y . 
2i]       1     D 


In  vector  notation 


6g         1      ,  ,    _J,_ 

6x       2 (Area)  iy23   y31   Y12  '  '  2A*    1 

(A. 6) 

where  A*  =  Area  of  triangle  1-2-3  and  B,  =  {y9o   y.,-,   Y-i?} 


In  a  similar  fashion 
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2A*   B2 


*7 


(A. 7) 


where  B2   =   {  x~2   x,_.   x~,  } 


and   x. .   =   x.  -  x . 
ID       ID 


Now  equation  (A. 5)  can  be  written  as 


6F 
6x 


2A*  Bl  t 


(A. 8) 


where 


i,      = 


l-2n v  n 


(4^-1  +  2n) 


(45,-D 


(l-4n)  v  ^1 


(453-  l  +  4n) 


4C 


d-2n) 


4C 


(l-2n) 


(l-4nP3 


d-4n) 


•U 


(l-4nP2   (i-4npl 
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In  a  similar  way 

*7  ' 


2A*  B2* 


(A. 9) 


Now  the  element  strain  vector  can  be  written  from 
equation  (A. 4)  as 


e  = 


> 
e 

X 

6u 
6x 

e 

y 

•  =   . 

6v 
6y 

•    = 

e 
.      xy   J 

6u        6v 
6y        5x 

f   •  u 
x   — 


f   •  V 

y 


f   •  u  +  f   «v 

I  y   -    x   - 


(A. 10) 


where 


6F  5F 

f   =  -j—  and  f   =  t— 
x    6x       y    oy 


e  = 


2A* 


Bn    0 


.  1 


B2  r 


B2    Bl 


±        0 

0    ij; 


•  < 


u 


v 


(A. 11) 


or 


e   =   b  •  > 


u 


v 


By  definition,  stiffness  matrix  of  the  element  is 


k  = 


vol  lime 


[b] 


[D]  •  [b]   d(vol) 


(A. 12) 


where  [D]  is  the  elasticity  matrix  and  k  is  the  element 
stiffness  matrix. 
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Ezuation    (A. 12)    can  be   expanded  as 


k   = 


T 
1 

0 

~»l 

0 

4 

Dn 

D12 

D13 

4A* 

, 

0 

0 

BT 
B2 

■d 

D21 

D22 

D23 

volume 

D 
31 

°32 

°3  3_ 

(1) 


(2) 


(3) 


~B1 

0 

• 

0 

B2 

_B2 

Bl_ 

(4) 


ty  0 


0  iJj 


d    (vol) 


(A. 13) 


(5) 


where  D. .  are  coefficients  of  elasticity  matrix,  as  shown 
in  (4.41). 

Multiplying  (2) ,  (3) ,  and  (4)  for  uniform  thickness 
of  the  element,  t 


k  = 


t 

• 

-  T      — 
V1         0 

• 

Cll  "•  C16 

4A*2   J 
are 

a 

0    V1 
(1) 

_C61  * " '  C66 
(2) 

>    0" 
0    ^ 

(3) 


d (area) 


(A. 14) 


where  C. .  are  the  coefficients  of  the  product  of  (2),  (3) 
and  (4)  forms  a  symetric  matrix  by  nature. 

Only  above  diagonal  coefficients  are  shown  here: 
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"1,1 


'1,2 


'1/3 


'1,4 


'1,5 


'1,6 


'2,2 


y23  ,(D11^23  +  D13X32 
y23  ,(Dlly31  +  D13X13 
y23  ' (Dlly12  +  D13X21 


y23  *  (D12X32  +  °13y23 


y23  *  (D12X13  +  °13y31 


=  y 


=  y 


23 


31 


'2,3  "  y31 


'2,4 


'2,5 


2,6 


31 


=  y 


=  y 


31 


31 


(D12X21  +  D13y12 

(Dlly31  +  D13X13 
(Dlly12  +  D13X21 


(D12X32  +  °13y23 


(D12X13  +  D13Y31 


(D12X21  +  °13y12 


C3,3  =  y12  -(Dlly12  +  D13X21 
C3,4  =  y12  *(D12X32  +  °13y23 


'3,5 


'3,6 


y12  ' (D12X13  +  D13y31 


=  y 


12 


(D12X21  +  °13y12 


C4,4  ~  X32  ' (D22X32  +  D23Y23 
C4,5  =  X32  ,(D22X13  +  D23y31 
C4,6  =  X32  ,(D22X21  +  D23y12 


+  x32  '(03^23  +  D33x32 
+  x32  • (D31y31  +  D33x13 
+  x32  •  (D31y12  +  D33x21 


+  x32  -(D32x32  +  D33y23 


+  x32  -(D32x13  +  D33y31 


+  x32  • (D32x21  +  D33y12 


+  x 


+  X 


+  X 


13 


13 


13 


+  x 


13 


+  x 


13 


*(D31y31  +  D33X13 
•(D31y12  +  D33X21 


(D32X32  +  D33y23 


(D32X13  +  °33y31 


(D32X21  +  D33y12 


+  x21  -(D31y12  +  D33x21 


+  x 


+  x 


21 


21 


+  x 


21 


(D32X32  +  D33y23 


(D32X13  +  D33y31 


(D32X21  +  D33y12 


+   y23    *(D32X32    +   D33y 
+    y..     • (D_0X,^    +    D-,y 


v    32    13  33jr31 

32X21    +    D33Y12 


+  y23  •  <D^xoi  +  d^y 
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C5,5  =  X13  -(D22X13  +  D23y31}  +  y31  * (D32X13  +  D33y31) 
C5,6  =  X13  *(D22X21  +  D23y12>  +  y31  *  (D32X21  +  D33y12} 


C6,6  =  X21  *(D22X21  +  D23y12)  +  y12  ' (D32X21  +  D33y12) 


Carrying  out  matrix  multiplication  in  (A. 14),  the  12x12 
coefficients  of  the  stiffness  matrix  can  be  determined. 
The  last  step  in  the  derivation  is  to  integrate  the  coef- 
ficients over  the  area  of  the  triangle,  multiply  by  thick- 
ness and  divide  by  four  times  the  area  squared.   The 
formula  given  in  equation  (3.16)  is  helpful  for  integration 
of  stiffness  coefficients.   As  the  coefficients  are  symmet- 
ric in  nature,  only  diagonal  and  below-diagonal  terms  of 
the  final  stiffness  coefficients  are  given  here  and  desig- 
nated as  Sk. .  where  i  and  j  refer  to  degrees  of  freedom. 

Note   A  =  A* 
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This  completes  evaluation  of  the  stiffness  matrix 
for  isoparametric  triangular  elements. 


Determination  of  Element  Strain  and  Stress 

In  general  finite  element  scheme,  stiffness  of  each 
element  is  joined  according  to  their  nodal  and  element  con- 
figuration to  form  the  structural  stiffness  matrix  which  is 
the  distribution  of  stiffness  coefficients  at  all  nodes. 
This  matrix  is  modified  for  given  boundary  conditions  and 
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the  system  of  equtions  is  solved  for  given  nodal  loads.  The 
solution  is  the  distribution  of  x  and  y  displacements  at  all 
nodes.  To  determine  strain  in  an  element,  nodal  displace- 
ments of  the  element  are  grouped  in  vector  u  and  v  for  x  and 
y  directions.  Then  equation  (A. 11)  can  be  used  to  find  ele- 
ment strain. 

To  evaluate  element  stresses,  only  the  elasticity 
matrix  [D]  and  strain  vector  e  are  required.   The  coeffi- 
cients of  [D]  matrix  are  the  same  as  used  in  equation  (A. 13). 
Element  stress  can  be  found  by  the  following  expression: 

{  a  }   =   [D]   •   {  e  }  (A. 15) 

If  incremental  analysis  is  used,  the  displacement 
vectors  u  and  v  ,  strain  vector  e  and  stress  vector  a  are 
increments  in  one  increment  of  load.   Total  displacements, 
strain  and  stress  can  be  obtained  by  algebraic  summation 
of  these  quantities  for  all  increments  in  the  analysis. 
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APPENDIX  -  B 
CURVED  BAR  ELEMENT 

Figure  B« 1  shows  a  segment  of  a  curved  bar  element 
1-2  of  uniform  thickness,  which  has  a  radius  R,  and  sub- 
tends an  angle  8  at  the  center.   E  is  Young's  modulus  of  the 
bar  material  and  I  is  the  moment  of  inertia  of  cross  section 
per  unit  length.   S  is  a  variable  point  along  the  arc  which 
makes  an  angle  4>  at  the  center  and  dS  is  the  small  increment 
in  S  corresponding  to  d<j> ,  the  increment  in  <t>  .   The  segment 
is  acted  upon  by  a  radial  force  Q,  tangential  force  P  and 
moment  M  at  each  end,  corresponding  deflections  are  v,  u and  6. 

The  following  sign  conventions  have  been  used,  v  is 
always  radial  and  positive  when  directed  towards  the  center, 
u  is  always  tangential  and  positive  in  the  counter  clockwise 
direction.   0  is  positive  in  counter  clockwise  direction. 

Total  moment  at  S  due  to  all  forces  and  moments  at  1 
is  given  by 

M  =  ML  +  PR  (1  -  Cos<j>  )  -  Q.R  Sine))  (B.l) 

1   f   2 
Strain  energy:   U  =  -^f       M  dS 

But  dS   =   Rd<j> 
R 


So  U  = 


2EI 


M2dc}>  (B.2) 

o 
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Center    of    Curvature 

1'U2 


%,\ 


Radius  =R 


FIGURE    B.I     CURVED     BAR    ELEMENT. 
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Displacement  in  any  direction  can  be  found  by 
Castigliano' s  theorem,  namely  -^=r  =  6. 

Or  •       1 

Displacements  at  node  1   are 


u. 


5U 
<5P. 


R 
EI 


M 


6M 
5P-, 


dcf) 


(B.3a) 


v 


5U_ 
6Q- 


R 
EI 


M 


6M 
6Q, 


d<j) 


(B.3b) 


6U_ 
5M, 


R 

EI 


M 
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(B.3c) 


Substituting  for  M  from  equation  (B.l)  in  above 
equations,  taking  the  derivative,  multiplying  and  integrat- 
ing, one  obtains  the  following: 
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c 

b 

a 
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R   2 
EI 
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QjR 
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(B.4) 


where 


a   =   6   -   Sintj) 


b   = 


Cos6  +  Sin_3  _  ! 


|3  _  2  sin3  +  Sin^ 


d  =  £■  - 


Sin  23 


e   =   Cos  3-1 
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Solving  (B.4)  for  P.,  Q,  and  M., 


M, 
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B     D 
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S.  R   I    R     L    R2 

p  R   3  R    6  R 
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B  =   - 
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C   =   ad   -  be 


D   = 


-  c 


E   =   ce   -  ab 


F  =  b   -  ad 


2         2, 
G   -   b(b  -  2  ^)   +   c(  |   -  d)  +  ^ 

Considering  equilibrium  of  the  segment  1-2  under  the 
system  of  forces,  the  following  set  of  equations  can  be 
written 

P.  +  P2  CosB  -  Q2  Sin3  =   0 

Q1   +   P2  SinB  +  Q2  CosB  =   0  (B.6) 

M,  +  P1R  (1  -  CosB)  -  Qj*   SinB  +  M2   =   0 

Solving  for  P„,  Q„  and  M»  in  terms  of  P, ,  Q,  and  M1 


Q. 


M. 


-  Cos  B       -  Sin  B     0 
Sin  B         Cos  B     0 

-  R(l  -  CosB)    R  SinB    -1 


rpll 

• 

Ql 

Ml 

I            J 

(B.7) 
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Substituting  in  (B.7)  for  P,,  Q..  and  M.  from  equation 
(B.5)  and  carrying  the  matrix  multiplication,  the  following 
is  obtained: 


M 
2 

where 


EI 


R3G 


[k21] 


u. 


V 


[k21]  = 


lll 


l21 


l12 


l22 


a31         a32 


and  a. .  are  defined  as  below: 


a.,  =  -  (A  CosB  +B  SinB  ) 
a12  =  -  (B  CosB  +D  SinB  ) 
a13   =   -  (C  CosB +E  SinB  )  •  | 


l13 


l23 


l33 


(B.8) 


L21 


l22 


23 


(A  SinB  -  B  CosB  ) 
(B  SinB  -  D  CosB  ) 

(C  SinB  -  E  CosB  )  •  i 

p 


l31 


32 


L33 


[A  (CosB  -  1)  +  B  SinB 
[B  (CosB  -  1)  +  D  SinB 
[C  (CosB  -  1)  +  E  SinB 


r 

6 

J 

•     R 

E 

e 

J 

•     R 

F 

] 

R 

•        

B 

write 
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Following  a  similar  procedure  for  node  2   one  can 


3    LK22 J 
R  G 


u. 


v. 


(B.9) 


where  [k0~]  is  given  by 


"22 


[k22]   = 


-B 


!» 


-  B 


B  R 


Ir 


£  R2 

6  R 


(B.10) 


and 


T"  [k2i] 

R^G 


u. 


v 


(B.ll) 


where  k~,  is  the  transpose  of  matrix  k,„. 

Combining  equations  (B.5),  (B.8),  (B.9)  and  (B.ll)  the 
stiffness  matrix  for  the  curved  bar  1-2  can  be  written  as: 
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M 

qi 

mi 

P2 

Q2 

M 

\                       J 

[    k  ] 


where 


[k]   = 


EI 


R3G 


or 


[k]  = 


SR(1,1) 


SR(6,1) 


(                \ 

ul 

vl 

el 

U2     " 

v2 

.    V 

kll    k12 


k21    k22 


6x6 


SR(1,2)  SR(1,6) 


SR(6,2)  SR(6,6) 


(B.12) 


SR(i,j)  are  the  stiffness  coefficients  in  the  local 
co-ordinate  system. 

The  final  step  is  to  transform  the  stiffness  matrix 
[k]  to  the  global  co-ordinate  system  x-y.   Figure  B.2  shows 
the  local  co-ordinate  system  u,v,  e   and  global  co-ordinates 
x-y.   Angle  <j>  is  always  measured  from  the  horizontal  line 
passing  through  the  center.   <j>  is  positive  counter  clockwise 
and  negative  clockwise. 

Expressions  for  transformation  u,  v,  9  and  x-y  system 
can  be  written  as : 
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FIGURE    B.2    TRANSFORMATION    OF   COORDINATE 

SYSTEM    FOR    CURVED   BAR    ELEMENT. 
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u      =      -   x   Sin<}>        +      y   Cos<J> 


v  =   -  x  Cost    -  Y  Sin  <J> 


=   6 


(B.13a) 
(B.13b) 
(B.13c) 


In  Matrix  form 


u 


v 


=   [  L  ] 


-     Sin<|> 

Coscj) 

0 

-     Cos<J> 

-   Sin<f> 

0 

0 

0 

1 

(B.14) 


where  [  L  ]   = 


Complete  transformation  matrix  for  element  1-2  can 
be  written  as 


'  ul 

vl 

el 

y 

U2 

V2 

62 

^ 

=   [L] 


*1 

91 
X2 
^2 


(B.15) 
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where 


[L]    = 


Sine)),         Cos<J>,         0  0  0 


Cose}),       -Sin*}).         0 


0 
0 
0 


0 

0 


0 
0 
0 


0      -    Sin<J>. 


0 
0 
Cost}). 


0      -   Coscf>2      -   Sin<J>2 


0 


Inverting  matrix   L   and  denoting  by  T, 

[T]   =   [L]"1 

pre-multiplying  [k]  matrix  by  transpose  of  [T]  and  post- 
multiplying  by  [T]  ,  the  stiffness  matrix  for  curved  element 
1-2  in  global  co-ordinate  system  can  be  written  as 


[SRG]     =  [T] 

6x6       6x6 


[k]  •   [T] 
6x6     6x6 


(B.16) 


where  [SRG]  is  the  stiffness  matrix  for  curved  bar  element 
in  global  co-ordinate  system. 
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APPENDIX  -  C 
INTERACTION  ELEMENT 

Figure  C-l  shows  an  interaction  element  with  four 
nodes   1  to  4  with  eight  degrees  of  freedom.   Initially  co- 
ordinates of  node  4  are  same  as  for  node  1  and  those  of  node 
3  are  the  same  as  of  node  2,  which  ensures  zero  thickness. 
The  local  co-ordinate  system  is  chosen  so  that  the  origin  is 
at  the  midpoint  of  the  element,  the  x  direction  is  parallel 
to  and  the  y  direction  is  perpendicular  to  the  length  of  the 
element.   The  stiffness  matrix  will  be  developed  using  the 
local  co-ordinate  system. 

If  the  interaction  element  is  acted  upon  by  a  vector 
of  force  per  unit  length  {P},  which  produces  a  relative  dis- 
placement vector  [w] ,  then  the  expression  for  total  potential 
energy  $,  of  the  element  can  be  written  as 


$   =•  =■ 


L/2 


-  L/2 


[w] 


{P} 


dx 


(C.l) 


where 


P  = 


[w] 

P 
s 

P 
n 


w 


top 


L.wn 


top 


w 


bottom 


w 


bottom 
n 


210 


8 


fe7 


df-i 


Y,cun 


Top 


6 

to 

©  —  5 


L/2 


Bottom 


•L/2 


44  -~  ^cuJ^JA  Zero 

&_3        3-  Th.ckness 


-I 


FIGURE  C.I     REPRESENTATION    OF    INTERACTION 
ELEMENT. 
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in  which  subscripts  n  and  s  refer  to  normal  and  tangential 
directions  respectively. 

The  vector  {p}  may  be  expressed  as  a  product  of  unit 
stiffness  and  displacements,  or 


(P)   =   [E] 


[w] 


(C.2) 


where  [E]  =  diagonal  interaction  stiffness  property  matrix 

E     0 
s 

0     E 
n 


in  which  E   and  E   are  the  stiffnesses  in  directions,  normal 
to  and  parallel  to  the  element. 

Substituting  equation  (C.2)  into  (C.l) 
L/2 


2 


[w] 


[E]    [w]   dx 


(C.3) 


-L/2 


The  next  step  is  to  write  the  relative  displacement 
vector  [w]  in  terms  of  nodal  displacements.   Along  the 
bottom  line  of  interaction  element 


bottom 


w 


bottom 
n 


w 


1 

2 


i-¥ 


1  +  2£    o 

Li 


i  -  2S   o 


i  +  22 

J-i 


' 

ul 

U2 

•     ■ 

U3 

U4 

I            J 

(C.4) 
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where  x  is  a  floating  point  along  x  -  direction  (Figure  C.l) 
Similarly  for  top  line  of  interaction  element 


w 


top 


w 


top 


n 


1 
2 


1  +  £* 


0 


i   ^x 

L 


1  +  ^ 


'   L 


.th 


f 

U5 

•  ■ 

U6 

■ 

U7 

U8 

,     J 

(C.5) 


where  u.  is  the  total  displacement  in  the  i    degree  of 
freedom. 

Now  the  relative  displacement  vector  [w]  can  be 
written  in  terms  of  equation  (C.4)  and  (C.5) 


[w]   =   J 


-A 


0    -B     0    B    0 


-A 


or   [w]   =   |  [D] 
where   A  =  1  -  ■= — 


0    -B 


{u} 


B 


and 


B  =  1  + 


2x 


ul 

U2 

U3 

« 

UM 

• 

U5 

U6 

U7 

U8 

(C.6) 


Substitution  of  [w]  into  equation  (C.3)  yields 
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L/2 


2 


|   (u}T  •  [D]T  • 


[E]  •  [D]  •  (u>  •  dx 


(C.7) 


-L/2 


Performing  matrix  multiplication  of  [D]     [E]    [D]  and 
integrating  over  -  L/2  to  L/2  yields  [k]  ,  which  on  sub- 
stitution into  equation  (C.7)  gives  the  expression  for  total 
potential  energy. 


2 


(u) 


[k] 


(u) 


(C.8) 


where  [k]  per  unit  length  (L  =  1)  is  given  by 


[E]  = 


2E 

s 

0 

Ec 
0 
-E 

C 

0 

-2E 

0 


0 

2E 

0 

E 

0 

-E 

0 

-2E 


n 


n 


n 


0 
2E 

c 

0 
-2E 

c 

0 

-E 


n 


0 
E 
0 

2E 

0 

-2E 

0 

-E 


-E 


n 


n 


n 


s 
0 


n 


-2E 


0 


s 
0   -2E 


n 


2E 


0 


s 
0    2E 


n 


E 


0 


s 
0     E 


n 


n 


-2E 

s 

0 

-E 

c 

0 

K 
o 

2E 

< 

0 


0 

-2E 

0 

-E 

0 

E 

0 

2E 


n 


n 


(C.9) 


n 


n 


8x8 


By  definition  of  minimum  potential  energy,  equation  (C.9) 
is  the  [8x8]  stiffness  matrix  per  unit  length  of  interac- 

tion  element. 

in  the  culvert  problem,  the  value  of  En  is  maintained 

very  high  as  long  as  the  element  is  in  compression  which  en- 
sures . essentially  zero  displacement  normal  to  the  pipe. 
Values  of  Es  are  also  very  high  before  a  limiting  value  of 
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shear  stress  to  normal  stress  (stress  ratio)  is  reached. 
For  values  of  the  stress  ratio  above  the  limiting  quantity, 
E   is  also  reduced  to  zero  permitting  slip  between  the  top 
and  bottom  faces  of  the  interaction  element. 

So  far  the  development  is  based  on  the  local  co- 
ordinate system  only,  depending  upon  the  element's  position, 
the  normal  and  tangential  directions  change.   The  next  step 
is  to  transform  the  stiffness  matrix  to  the  global  co- 
ordinate system,  which  follows  a  standard  procedure  of  co- 
ordinate transformation  as  described  in  Appendix  -  B,  equa- 
tions similar  (B.13  to  B.16). 
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APPENDIX  -  D 

Toct 
DERIVATION  OF   

°oct    failure 

In  equation  (5.2)  ratio  of  octahedral  shear  stress 
to  octahedral  normal  stress  at  failure  is  presented.   The 
derivation  of  equation  (5.2)  follows. 

c  ,  a2  and  c?3  are  the  principal  stresses  at  failure 
for  a  soil  whose  angle  of  friction  is  tj>. 

a,  +■  a~  +a-, 
12    3  (D.l) 

aoct  -       3 


*oct-  lAl-^   +  (02-°3)   +  (a3"al} 


(D.2) 


For  plane  strain  conditions  from  the  generalized  Hook's  law, 
a     may  be  represented  as 

c2  =  *(a1  +  a3)  (D-3) 


where 

tjj   =   v 
Substituting   for   a2   from  equation    (D.3)    into    (D.l) 

and    (D.2) 

(°1    +    °3)(1    +    *j  (D.4) 

°oct  3 

and 


2    H>2   -  ib  +  1)    -   6   0,0,  (D.5) 


T  =   \    /2(o,    +   a7)2     (iT   -   *   +   1)    "    6    °ia3 

oct  3    /  -L  J 


Dividing   (D.5)  by  (D.4) 
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Toct    /2((J1  +  a3)2  ^   ~   *  +  1)  "  6ola- 


oct 


[a1   +a3)  (1  +\p) 


(D.6) 


^4.^-1       1    1  +  Sine}) 

At  failure  — —  =  -= =-^ — f  =  N^ 

a    1  -  Sintj)    4) 


(D.7) 


From  equation  (D.7)  one  can  write  the  following 


(a1  +  o3)  =  (N   +  l)-o3 


al   =   %a3 


al  a3  =   Vb' 


Substituting  for  a..  ,  a,  etc.  from  above  identities  into 
(D.6)  and  simplifying,  one  can  write 


oct 
7oct 


i  +i> 


failure 


6N 


2(ijj   -  ij;  +  1)  - 


(N.+  1) 


(D.8) 
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APPENDIX  -  E 

DERIVATION  OF  TANGENT  MODULUS  (Et)  AND  TANGENT 
POISSON'S  RATIO  (vfc) 

For  a  three  dimensional  state  of  stress 
a,   =  major  principal  stress 

a        =   intermediate  principal  stress 

o^   =  minor  principal  stress 

e    =  major  principal  strain 

e    =  intermediate  principal  strain 

e    =  minor  principal  strain 

E    =  Young's  modulus  of  the  material 

v    =  Poisson's  ratio 
Generalized  Hook's  law: 
E'el  =  °1  "  V°2  "  Va3 
E'£2  =  °2  "  V°l  "  V°3 

E*£3  =  °3  "  V°l  "  V<J2  (E.3) 

For  plane  strain  condition  e2  -  0  and  from  equation 

(E.2) 


(E.l) 
(E.2) 


(E.4) 

2  "    VKU1    'u3; 


o^  =   v(on  +o,) 
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Substituting  a-  from  equation  (E.4)  into  (E.l)  and 
(E.3)  and  simplifying 


1  -v' 


a,  - 


1    1  -  v  u3 


(E.5) 


1  -vz  r 


E 


a->  - 


3     1  -v   ul 


(E.6) 


:-,         a,  -  U )  a -j 

-  K  - 


a,  -  h )  a, 

3     1-v    1 


(E.7) 


or 


1  -v 


al  -Ka3 

a3  "  Kal 


or 


1  -  Ka. 


v  = 


(1  -  K)  (a1   +a3) 


(E.8) 


Substituting  for  K  from  (E.7)  in  (E.8)  and  simplifying 


v  = 


Q1     £3  -  £l  O3 

(e3  -  e,)  (oT  +  oT) 


(E.9) 


To  get  tangent  Poisson's  ratio,  v.  incremental 
values  of  stresses  and  strains  have  to  be  used.   If  A  denotes 
incremental  values,  the  expression  for  tangent  Poisson's 
ratio  is  given 


Aa,Ae3  _  Ae-,Aa-, 
vt  "   (Ae3  -  As.)  (Aa,  +  Aa3> 

Subtracting  (E.6)  from  (E.5)  and  simplifying 


(E.10) 


or 


Substituting  for  v  from  equation  (E.8) 
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E  = 


1  -  v 
£1  "  e3 


ai  "  a3 


(E.ll) 


or 


E  =  (1  +v) 


al  ~a3 
£1  "£3 


(E.12) 


E  =   U  + 


°1  "  K  °3 


(1  -  K)  (a1  +  a3) 


°1   ~    °3 
Cl    ~    e3 


(E.13 


a1    (2  -  K)  +  a   (1  -  2  K) 
(1  -  K)   W~i~a^) 


al  ~  °3 
el  "  £3 


(E.14) 


E  = 


al  -  P3 
e  -  c3 


°1    £  3  ~  £1  03 
(a1  +  a3)  feT  -  e_l 


(E.15) 


To  obtain  tangent  modulus  E ,  ,  incremental  values  of 
stress  and  strains  are  used  in  equation  (E.15) 


(Act,  -  Aa3) 
Et  =  (Ae,  -  Ae0) 


1  - 


Aa, Ae_ 


Ae,  Aa_ 


(Aa,  +  Aa,) (Ae,  -  Ae_) 


(E.16) 
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